{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "32ce4267",
   "metadata": {},
   "source": [
    "# 子空间方法的频率估计"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10024526",
   "metadata": {},
   "source": [
    "创建一个长度为 24 个样本的复值信号。该信号由频率为 0.4 Hz 和 0.425 Hz 的两个复指数（正弦波）和加性复高斯白噪声组成。噪声的均值和方差为零0.2的平方在复数白噪声中，实部和虚部的方差都等于总方差的一半。\n",
    "因为添加有随机高斯白噪声，每次运行结果不一样。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c6c155a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABYoUlEQVR4nO2dd5wcZ3n4v8/W6/0k3amdZEuWZVlukgu2sVzoprdA6AkmQAKEAAYMCTWh5JcQIAnNodlgerMB4yb3IsmWJcuyeruTTrretz+/P2Zmb6/s3u7etju9389nP7c7Mzvzzt7MPO/TRVUxGAwGgyERV7EHYDAYDIbSwwgHg8FgMEzBCAeDwWAwTMEIB4PBYDBMwQgHg8FgMEzBCAeDwWAwTMEIB4PBYDBMwQgHQ0kgIodFZExEhhNerUUczz0ioiLiSVj2eRHZKSIREflMGvu4UEQesM/lpIh8MGHd+SLyoIgMiEi7iHx60nffICK7RWRIRJ4VkVclrBMR+YKIdNjf3ywi58ziXFVEzpy07DMicku2+zTMfYxwMJQSL1fVqoTXcWdF4kM634jIXwPeaVbtBz4G3JHGPpqAPwPfBhqBM4G/JGzyE+ABoAG4CnifiLzC/u5i4Bbgw0AN8FHgJyKywP7u64F3AVfa338U+HFGJ2kwzIARDoaSxZ7Rvl9E9gH77GXXi8h2EekXkUdEZH3C9jfas+khEdkjItdmccxa4F+whMAEVPWHqvonYCiNXX0YuFNVb1XVoKoOqeruhPVtwK2qGlXVA8BDgDP7XwL0q+qf1OIOYAQ4w16/AnhIVQ+qahRLkKxNOIfNtmbxiK21/EFEGkXkVhEZFJEtItKWwW/ysUkaXVhEfpDu9w1zEyMcDKXOq4BLgLUicgHwf8B7sGbj3wZ+LyJ+ETkL+Htgo6pWAy8CDgOIyJttYZLstSzheP8K/C/QOctxXwr02g/oU/YDOvE4XwPeJiJee+yXAXfb67YCu0XkFSLitk1KQWCHvf424AwRWS0iXuDtWFpKIn8FvBVYjCVUHgW+j6Vp7MYSgGmhql9xtDngbKAL+Fm63zfMUVTVvMyr6C+sB/kw0G+/fgsocE3CNv8LfH7S9/ZgmWXOBE4B1wHeLMewAdgOeLBm9gp4ptnuFuAzM+xrr30eG4Ey4OvAwwnrn4dlporYx/nspO//jf17RIBR4GUJ63zAf9nfiwCHgBUJ6zcDNyV8/n/AnxI+vxzYnvBZgcGE374fCAC3TBpTObANuLHY14t55f9lNAdDKfEqVa2zX6+ylx1LWL8c+KfEWT+wFGhV1f3Ah4DPAKdE5LZMHNoi4gL+B/igqkZmfyqMAb9R1S2qGgA+CzxPRGpFpAFrpv85LMGxFHiRiLzPHst1wFeATViC4CrgeyJyvr3vf8YSOkvt738WuFdEKhKOf3LSWCZ/rpo03gsTfvs64EvTnNPNwB5V/XLav4JhzmKEg6HUSSwbfAz4YuJDTFUrVPWnAKr6E1W9AkuIKPBlsBzMk2zmk1/LsBy/G4CfiUgnsMU+ZruIXJnFuHdMGnvi+5VAVFV/pKoRVW3HMhW91F5/PvCAqm5V1ZiqbgEex9KKnPU/U9V2+/s/AOpJ8DvkGhH5OLAaS6MxnAYY4WCYS3wX+DsRucQO56wUkZeJSLWInCUi14iIH8skMgbEANRyCleleB0FBoBWrAfv+Yw/qC/CejBj+wfKsO4bj4iUiYg7yVi/D7zaDln1Ap/GciIPYJmcxPaFuERkEfBGxn0KW4ArHU3B9rVcOWn960Vkof39t2JFV+2fzY+bDBF5CfAB4NWqOpaPYxhKDyMcDHMGVd0KvBv4JtCH9TB8h73aj2UK6cZyJi8APpHBvlVVO50XltMV4KSqhuz338USOm8CbrLfvxVARK4UkeGE/d0LfBIr7PUUlk/kzfa6QeA1wD/a57EdeAb4gr3+fizz2C9FZAj4FfCvquqEwn4ZeNr+Xr+9n9eqan+655shbwSasZzkjrb1rTwdy1AiiKpp9mMwGAyGiRjNwWAwGAxTMMLBYDAYDFMwwsFgMBgMUzDCwWAwGAxTKFgxs3zS1NSkbW1txR5GxoyMjFBZWVnsYRQUc86nB+ac5wbbtm3rVtXm6dbNC+HQ1tbG1q1biz2MjNm8eTObNm0q9jAKijnn0wNzznMDETmSbJ0xKxkMBoNhCkY4GAwGg2EKRjgYDAaDYQpGOBgMBoNhCkY4GAwGg2EKRjgYDAaDYQpGOBgMBoNhCvMiz8FgMBhKlUPdI/xi6zEUeN1FSzijeXITvtLECAeDwWDIE5v3nOKGH28jFlNE4PsPH+J7b9vIFauaij20GTFmJYPBYMgD24/1c8OPt7FqQRWPfOIaHv74NbQ1VvK+W7dxaihQ7OHNiBEOBoPBkGPGQlE+/LPtNFf5ufVvL2FBdRkLqsv4n7++kEA4xv+7c2+xhzgjRjgYDAZDjvnf+w9wsHuEr7xuPXUVvvjylc1VvPmSZfzqyXZODJR2O24jHAwGgyGH9I6EuPnBg7zs3BYuP3Oqb+Gdl7cRiSm/fep4EUaXPkY4GAwGQw759gMHGA1H+dB1q6Zdv7yxkouW1/Obp9pR1QKPLn2McDAYDIYccWoowA8fOcyrzl/MqoXVSbd71QWL2XtymGdPDBZwdJlhhIPBYDDkiP/dfIBwVPngtdNrDQ4vO7cFEbhn96kCjSxzjHAwGAyGHNAbiHHrY0d53YVLaGtK3RGuodLHOa01PLy/u0CjyxwjHAwGgyEH/OFAGEX5h2vPTGv7553RxFNH+xkLRfM8suwwwsFgMBhmyZGeER5oj/DGjUtZUl+R1need0YjoWiMrUd68zy67DDCwWAwGGbJV+7cg9sFH7gmta8hkY1tDXhcwiMHevI4suwxwsFgMBhmwfZj/dyx4wQvafOyoKYs7e9V+j2saalmZ/tAHkeXPUY4GAwGQ5bEYsoXbn+WpiofL17hzfj761preeb4QEnmOxjhYDAYDFly6+NH2Hqkj4+9aA3lHsn4++sW19I/Gqajv/RKaRRFOIjIV0XkORHZISK/EZG6hHXrReRREdklIjtFJH09zWAwGApEe98oX/rTc1y5qonXb1iS1T7WLa4F4JmO0kuGK5bmcBewTlXXA3uBTwCIiAe4Bfg7VT0H2ASEizRGg8FgmBZV5ZO/eQYF/vXV5yKSudYAsGZRNW6XsOt46fkdiiIcVPUvqhqxPz4GOGL3hcAOVX3a3q5HVUszCNhgMJy2/G77cR7Y28XHXnQWSxvSC12djjKvm1ULqnimo/SEgxTbESIifwB+pqq3iMiHgIuABUAzcJuqfiXJ924AbgBYuHDhRbfddluBRpw7hoeHqaqaGy0Dc4U559OD+XzOwyHlEw+N0lzu4lOXluGytYZsz/lbTwc40B/jq1dlL2Sy5eqrr96mqhumW5e3NqEicjewaJpVN6nq7+xtbgIiwK0J47kC2AiMAveIyDZVvWfyTlT1O8B3ADZs2KCbNm3K+Tnkm82bNzMXxz0bzDmfHsznc/7Er3cwGhnjm2+/nLWtNfHl2Z7zU+G9PH7vPi674kr8HncORzo78iYcVPW6VOtF5B3A9cC1Oq6+tAMPqGq3vc0fgQuBKcLBYDAYCs2+k0PctuUY77p8xQTBMBtWNleiCkd7RlNWci00xYpWejHwMeAVqjqasOpO4FwRqbCd01cBzxZjjAaDwTCZb963n3Kvm/dfnV79pHRY2WSZog50jeRsn7kgb5rDDHwT8AN32V7+x1T171S1T0T+A9gCKPBHVb2jSGM0GAyGOMd6R/nD08d595Uraaj0zfyFNGlrsnwNh7qNcEBVk4pdVb0FK5zVYDAYSoZfbD0GwDsub8vpfqvLvDRX+znUPZzT/c4WkyFtMBgMMxCNKb/c1s6Vq5ppqS3P+f5XNlVysMTMSkY4GAwGwww8cqCb4wMB3rBhaV72v7K5suTMSkY4GAwGwwzcuauTCp+ba89ekJf9L6mvoGckVFKNf4xwMBgMhhSoKvfuPsUVZzZR5s1PHkJrnVVC7vhA6RTgM8LBYDAYUvBc5xDHBwJ50xqAuB/jRH8gb8fIFCMcDAaDIQX3PncKgKvPyp9wWFxnCQejORgMBsMc4bGDPaxZVJ1Rl7dMWVhThggcL6G+DkY4GAwGQxIi0RhPHuljY1tDXo/j87hoqvIbs5LBYDDMBZ49MchIKMrGFfkVDgCttWXGrGQwGAxzgS2H+wDY2Faf92O11pUbs5LBYDDMBbYe7mVJfXlesqIn01JbzomBAMXuseNghIPBYDAkYUf7AOcvrSvIsVrryhgNRRkci8y8cQEwwsFgMBimoXckREf/GOcuri3I8RbVWtFQnYOl4ZQ2wsFgMBimYddxq6/zugIJh8ZKPwA9w8GCHG8mjHAwGAyGadh1fBCAc3LU8W0mmqqsHhHdI6GCHG8mjHAwGAyGaXimY4Al9eXUVeSusU8qGqsszaHXaA4Gg8FQuuw6Psi61sKYlADqyr24BHqM5mAwGAylSSAc5XDPCGtaqgt2TJdLaKj00z1shIPBYDCUJAe7RlCFMxdUFfS4TVU+45A2GAyZs7N9gP+8ay/37+0q9lDmNQe6rH7OZzQXVjg0VvlKxqzkKfYADAZDetyz+yQ3/Hgb0ZiVQfvRF53F+68+s8ijmp/sPzWMCKxoqizocRsr/exo709r29FQhDKPG5dL8jIWozkYDHOAgdEwH/3lDs5uqWbbp67jFee18u9/2cO2I33FHtq85EDXMEvrK/LW+S0ZjVU+etLwOXzz3n2c+5m/8IL/vJ8TeSrWd9oLh+P9Y/GZmMFQqvzgkcP0joT48mvX01jl599ecy6NlX6+dvfeYg9tXrL/1HDB/Q0AjZU+hoIRAuHkvaS3Henj3/+ylw3L6zk5GOTjv9qZl7Gc1sLh0QM9XP7le3lof3exh2IwJCUSjXHL40e4+qxmzrFDKyv9Hv72yhU8uK87bTOEIT2iMeVQ9whnNBfWpAQJuQ4p/A7fe/AgteVevv/OjXzv7Rv419ecm5exnNbC4cLlddSUefnVtvZiD8VgSMrDB3roGgryxo3LJiz/60uWUVPm4eaHDhVpZPOTjr4xgpFY0TQHIKlpaSQY4Z7dp3jNhYup8Hm4dGVjvMVorjmthYPf4+Zl61u469mTBCPJ1TiDoZjc9WwnFT43V69pnrC8uszLy89r5c5dnQwHS6OS53ygWJFKMK45dI9MH8764L4uQtEYL1y7KO9jmVE4iMgCEXm1iLxfRN4lIheLyLwRKtedvYCxcJQnDvUWeygGw7Q8tK+bS1c24vdMdY6+5sLFBMIx/vxMZxFGNj/Zf6p4wqHB1hz6kpiVHjvYS7nXzYYCNB9K+pAXkatF5E7gDuAlQAuwFvgUsFNEPisihalIlUcuW9mEz+3ioX3G72AoPY71jnK4Z5Qrzmyadv2Fy+pZ3ljBb5/qyNsY7t/bxY2/3MGjB3rydoxS4kDXMI2VPuorC1NTKZGaMiu7YCgwvSb41NE+1i+pxevO//w8VZ7DS4F3q+rRyStExANcD7wA+FWexlYQyn1u1i2uYasJCTSUIA/bwRJXrppeOIgIrzivlf++bz/dw0GabLNErnjyaB/v+sEWojHlN9s7+PV7n1ewEtbFYv+pYc4ogr8BLFMhwOBYeMq6YCTKruOD/O2VKwsylqTiR1U/Op1gsNdFVPW3qjqnBYPDxrYGdrYPpAwfMxiKwaMHe1hQ7U/pHH3Z+hZiSl5MS//vL3torPTxwEevpsrv4f/9ZU/Oj1FqHOgaLopJCcDncVHudTMYmCocDnaNEIkpawtUQjylbiIiV4nIevv9G0TkmyLyjyKS2+lJkTlvaR2haIx9J4eLPRSDYQI72we4YFkdIsmzYM9aWM0ZzZXcseNETo99pGeEh/f38PbntbGssYK3X9bGfXu6ONY7mtPjlBI9w0H6RsNFiVRyqC7zTNsqdO/JIcD6fxeCVD6H/wa+AHxPRG4B3gw8A1wI/F9BRlcgnGYeTuenTPjtUx28/BsP8anf7jSahyGnDAXCHOwembFstIjwsvWtPH6oh1NDuWsxeecuSxN5xXmtgOX8BrhjZ26FUC4IRqI5SWYdd0YXPsfBoabcy1Bwquaw9+QQHpcUrKRHKs3halW9Eng+lkP6tar6LeBtwPpCDK5QLK2voNrviXd+SpdHD/TwoZ9tp280xC2PHeVztz+bpxEaTkeeta/HdUtmtvFfb5uW/rQzd6al+/d2sWZRNUsbKgBY2lDB2pYaNu85lbNj5ILbdxznvM/+hed/5b74wz1bDnSNAIWvxppITRLN4WDXCMsaKvB5ChMsmuooAQBVDQBHVDVqf1Zgqlibw7hcwlmLqtljq23poKr82592s6S+nLv+8SredfkKfvrEUfafSn8fhtJkNBThC7c/y1tvfpzfbc9fFNBM7Oywexin0XBm9cJq1iyq5pc5SugMR2M8dbSfS1c2Tlh+xaomnjzSz1ioNLTk3pEQH/3FDpbWVzAcjPCRXzyN9YjKjv2nhin3ummtzU9iWTrUlHun9Tkc6xuNC+pCkEo4LBCRD4vIPyW8dz43p/jenOTMBVUcyGDWsaN9gB3tA7zn+Ssp97n5+2vOxOtycevj0/rwDXOEWEy54UfbuPnhQxzsGuGDt23nF1uPFWUsz3QMsKimjObq9Fx8b7p4GTs7BqYU4wtFYhk/MHefGGQ0FJ0ST3/ZykZC0RhPHSuN6L5bHzvCWDjK/77lQj7+kjVsP9bPowezD7k90DXMyubKvFU6TYfqMu+00UpHe0ZZViLC4btANVCV8N75/L38D62wnNFcRc9IKGnyyWTu2HkCr1t45QWWHbah0scL1i7kd9uPm0J+c5jbthzjof3d/Ourz+X+j27i0pUNfO72Z1PWuskXu44Psm5x+pEpr9+whIZKH1/+83NEY8rJwQAfvO0pzv7nP3PFl+/jyaPpP9Cf6bBMWuctqZuw/Pyl1ucd7Zn75/LBHTtPsLGtnjMXVPPqCxZTU+bh109mr+3tP1W8SCWHmjLPlDyHgdEwg4FIaQgHVf1sqlfBRlggHBvj/q70tIe7nj3JZWc0UWPHJQO8eN0iekdCbC+RWdVc4ZH93Vz11fu44sv3ct9zxbNnh6MxvnHvPi5aXs9fbVyKx+3ic69cx1Agwo8ePVzQsUSiMQ73jHDmgvQjUyp8Hj7+4jU8caiXV3zzIa7598386ZlO3nTxUlwueM+PtzEwzYx0OnafGKTa72FJ/UTzSn2lj+WNFWw/2p/J6eSFjv4xnuscipeSKPO6ue7shdz17EnC0VjG+xsORujoH2P1wiILB9uslKjtHeuzIsQm/z/ySapopa+nehVshAWizY4AONIzc5je8f4xDnWP8PxJiUnPX92M2yVs3mO6dKXLiYExbvjxNgSo8Ln5u1u28VxnZoEBueLuZ09yYiDAe686Ix46unphNZvOauYnjx8tqEbY0T9GOKqszDBq5vUblvCZl68F4AVrF3LXPz6fL7zqXP77zRfSNRTk1sePpLWf3ScGWdNSPW0I7brWWnYX6X+UyBa75M0VCffhtWcvZGAsnHFwCcA+2+e4ukChosmoKfMSjiqB8LiA6xywotBa8lRkbzpSmZW22a8yrPDVffbrfGBWeeUi8nkR2SEi20XkLyLSmrBuk718l4jcP5vjZMLiunJE4GgaMdxbDlsX5WRnXW25l3WtNaZOUwZ84979hCIxfvSuS/jJuy+l3Ofmi3fsLspYfv1UBwtr/Fy9ZsGE5a+7aAmnhoIF/b8etKNmVmYYtigivOPyFdzxgSv52l9dwPJG6/vrl9Rx2cpGfvrE0bT8D/tODSd9SJ6xoIpjvaNFD93edqSPKr9nwjg32j6SrYcz/1/t6bSEw5pFxa0KVG2X0Eh0Sp+0Q5QX1ZQVbBypzEo/VNUfYoWtblLVb6jqN4BrsQTEbPiqqq5X1fOB24F/BhCROuB/gFeo6jnA62d5nLTxeVy01pZztGdkxm2fPNJHhc/N2S1TL6INbQ1sP9ZvqrymwXAwwm+f6uBVF7SyrLGCpio/773qDB7c151VzslsGA1FeGBvFy8+ZxHuSc7Ia9csxO9xcdezJws2Hqcy6Moc2r9fcX4rx3rH2DtDsmffSIiBsXDSePozF1QRUzjUPfO9kk+ebu/nvKW1E/5fC2rKWN5YEZ/AZcKek0OUe90FNd1MR0351BIaJwcCiEBTVeHqPaUTMFsPJD4Fq+xlWaOqiTpfJeBMZd4M/Nop26GqBTVAL2uoSEtzeOb4IOe01kx5iABctLyeYCQWn4WkQyymfPeBg7zhW4/ylT8/V/QZWaG4+9mTjIaivGHD0viyv7p4GWVeFz99orBRXw/s7SIYifGic6aWQi73WVUwHzlQuOKMB7tHqKvwxqt05oKrz7I0ovtmyFM4aD/0kwoHW2DNNqcgkW1Hennpfz3IR37xdFr+gmhM2dM5xNnTzPI3LG9g6+G+jCO09p4cYvXCqqJGKsF48b3BBKf0yUGrbpanAAX3HNI50peAp0TkByLyQ+BJ4F9ne2AR+aKIHAP+GltzAFYD9SKyWUS2icjbZnucTFjaUM6xvtT9WKMx5dnjg/GOXJNZs8hScZ/LQDh84979fPGPu+kdDfE/mw/wvlufJHYaRDzdt+cUTVU+Llw2PteoLfdyzZoF3LnrZEF/g7/sOkldhZeLVzRMu/55ZzTxXOcQ3cPT19nPNQe7hjM2Kc3EotoyVjRVzth3+rAtHNqSHH9lcyUiuRMOgXCU997yJM91DvLLbe1pNS862jtKMBJj9aKppq+NbfX0jIQy1mz2dA5x1jT7KzRxzSHBrNQ5GCioSQlSV2UFQFW/LyJ/Ai6xF92oqjOmYYrI3cB0HSluUtXfqepNwE0i8gng74F/scdzEZbpqhx4VEQeU9UpjXJF5AbgBoCFCxeyefPmmYY0I8G+EN1DYe6+9z48SWYPnSMxxsJR3IPH2bx5quM5porPBXdveZYFwwdSHm94eJhf/PFevv7gGJe2uHnP+hj3HPVxy+5TfP4nd3PVEm/K789FhoeH2bx5M6rKvbtGObfZwwMPTHQtLXNF+ONQkO/97l5W1+e/wbuqcu+zY6xpcPHQgw9Mu035gKXN3fyHB7mkZcbbZgLOOWfCcx2jrGty5+S6TqTVF+Tx/aPcd999Ses1PbAvhACHdm7hWJL7oKlMeHTXQS7wHp92fSbnfH97mFNDIW7cWMbtB0N86949rI4dxZWintS2k9aseqRjL5sn3WeRIUvzuO2ux3hea3r/q8Gg0j0cwjV0KuvfPJv/83QcH7bG//i2HcgJa/wHT4zRUCY5vx5SoqrTvoC2ZOvs9QIsSbVNOi9gGfCM/f7jwGcT1t0MvH6mfVx00UWaC376+BFdfuPterRnJOk29z53UpffeLtuOdSTdJtXfONBffN3H53xePfdd59+9ve79MxP3qHH+0dVVTUWi+mr/vsh3fiFu3QkGM78JEqc++67T1VV950c0uU33q63PXFkyjZDgbCuuumP+tnf7yrImI72jOjyG2/XHz16OOk24UhU1/3zn/Xjv9qR8f6dc06XkWBYl994u37z3n0ZH2smfvTIIV1+4+3a3jeadJuP/Hy7XvzFu1Lu553ff0Jf9J/3J12fyTm/9ebHddNX79NYLKa/fapdl994u2470pvyO/9z335dfuPtOhSYeo+EIlFd9ck/6r/e8WzaY7hrV6cuv/F2ffxg8vt6JjL9Pyfj5MDYlOvx4i/epR/9xfac7D8RYKsmea6mMit9VUR+JSJvE5Fz7I5wy0TkGhH5PPAwcHY2AklEViV8fCXwnP3+d8AVIuIRkQosbaVgoStOmFjnYPLiZUdsVdWJApmONYtq2H1iaEabZySm/HZ7By9Yu5AWO11fRPjUy87m1FCQWx+bv9nW245YDsOLlk8141T5PVx+RmPBavg4zsuNKbpredwuzl9Wx472/ryPp8M2bebDMbrGDqLYm6JUzPGBMVpnCJk8o7mSg90jszb9jQQjPHqgmxesXYiI8PxVzbgE7p8hHPxIzwhNVT6q/FM1A6/bxaqFVezOwLS77WgfHpewPo06VvnGMSsN2WYlVaV3JFTw5kOpopVeD3waOAv4b+BB4PfAu4E9wDWqeleWx/2SiDwjIjuAFwIftI+5G/gzsAN4Avieqj6T5TEypqXWsukd70/udzjcM0qlz50yauCsRdX0joTomsE+va8vRu9IKF710uGi5Q1sbKvnlsePzFvfw86OAar9nqTVLy8/s4mD3SPx+O58suVwH9VlHlbPkHC2trWGvSeHCEUyT7DKhPb+/AkH5xz3pRIO/YEZhcOS+gpCkVjSXsfp8tTRfsJR5XK70119pY+zW2pmzOY+MkMpiTWLanjuRPq5Do8d7OGcxbWUefNvxpwJv8eF1y3x4ntDwQjhqNJYKsIBQFWfVdWbVHWTqp6lquer6ptU9Ra1CvJlhaq+VlXXqRXO+nJV7UhY91VVXWuv/1q2x8gGRzikeiAd7hmhrakyZX39NS22U/pE6pnL011RfG4XV6yaWqrqLZcu50jPKA/tL+32pU6k1T/+bHvKB85k9p4cZvWi6ZOsYDyH5NGD+T//rYd72bC8fsYolXNaawlHlX15Lq7oTE5mekBnQ22FlwXV/qThrKpKR/8Yi2c4tjO2E/2zE95bj/QiAhcsq4svW7+kjh3tAyk176O9oym197Nbqjk1FKQnjQCC7uEg24/1c/VZpVEyTkSo9HsYCVrCoXfYKt3SUFnYNjozRiuJiFtEXiEiH0govvfhQgyu0FT5PZR73XQNJb+gjvSM0pbiooTxUL/DM+RMPNcb5cLlddOqxi9et4iGSh8/K1LRt3T5/iOH+eIfd/Obpzp4681PMBycvvdtIqrKPjtsMBlrW2qoLffyyP789i3uGwmx79QwG9qmj1JKZLzvR36zgzv6xvC4hAXV+YlOWdlcGY9ImkzvSIhQJDZjZExr3cxadjo8fayf1QuqJ5ShOXdxLQNjYY71Tr/vSDTGiYExlqbQrJwcpHSiBm9/+jiqVkZ5qVDp8zASsu6lHruuV0lpDjZ/AN4BNDJefK/48V55QERorvYnNQdFojGO9Y7S1pS6+FVztR+/x8XRFKU4hoMRjgzGuDjJQ8nvcXP9+hbu2X0yPoMoNcZCUb5+zz6uWt3Mr957GZ2DAW5LIz+hy+62lapMgcslXLqygUfy3NTeCevcmIZwaGuspMLnjvdZyBcd/WO01JVNm0eTC1Ll8zjX/sKZhIPtIzs+S7Pfc51DU9peOpp3Mr/IqaEgMYVFKcpqOyHlu2cwLfUMB/nug4e4YFld0vD0YlDpd8fve6cYaC5zXtIhHeGwRFVfo6r/ovO48J5Dc7U/aSx7R/8YkZimVGfBEjJLGyrixbKmY8exfhS4KMVD6fr1rQTCMe7eXbjM3Ez4484TDIyFee+mM7hoeQMbltfzi60z9xNw2rHOVMPm0pWNdPSP0THL2WkqthzpxetOzxHpdglrFlXnXzj0jeW1n8DS+gpODQWn7cngaM0zlQmvq/BS7nXPSnPoHw1xYiAQf5A7nNFkaZQHu6c3fZ1w6gzVJhdgjVV+mqr8KR3vI8EI7/rBFnpGgnz6+rWZDj+vWGYl6//TW8LC4U8i8sK8j6REaK7yJzUrHbY1gZnMSuDMzpLfOE4kxdppSnA4bFhez6KaMv7wdOm1ZQSrjWRLbRmX2IljLzm3hT0nh5KaLBycG3bVDNUvndn8lgxqGjlJiicG0ntobT3cx7kZOCLXtNSwtwA+h8V5LOGwrNHSfNunmbykKxxEhNa6slkJB6eKwOTEs9oKL42Vvnh9qck4PsFFKYQDwOqFVexJ4lsJR2O879Yn2dkxwDfedOGERMxSoMo/blZyKunWVRQ27ykd4fAY8BsRGRORQREZEpHil2TME83VyYVD54DjKJzZFrysoYL23tGkTrU9nYPU+FLfhC6X8OJ1i3hwX1fJdN5yCEViPLivm2vPXhB3Kl9jF6x7eIYyE3tPDlFf4aW5KvUD6OyWGqr9Hp5Is07Oc52DvOA/7+elX3+Qy/7tXj75m50pSzEEwlF2tPenZVJyWNFYSf9omP7R/PR3CEdjdA4GWJLH6ptL6i3hMJ1mm65wAMspPRuz0hHbtLWyaeokYWVzZVLh4Aj+VJoDWJrpvpND00b8ff2efdy/t4svvvrckvI1OFT4xs1Kg4EwLrH8EIUkHeHwH8BlQIWq1qhqtaoWt2xhHmmu9tM3Gp42XLFzwLpx0nEULqkvZygYoX90+vr5z3UOsaR65p//BWsXEozESi5qaWfHAGPhKFecOV4uuc0unrft8EzlGUZZ2VyVMuILLDPOhcvr09IcTg4GeNvNTzASjPDV163nXZev4CePH+Uzv9+V9Ds72gcIRzUtZ7SDU1IiX0XnOgcCxDQ/kUoO41F5UydBXUNByr1uKn0za1KtteWz0hzae0dxCbRMM9la2lAxrWYD1v/a73FRW556Jn3WompGQ9EpZsljvaN86/4DvPqCxbzp4mVZjz+fJJqVhgIRqsu8Ba/5lI5wOIaVwTw/A+4n4dj1ppsZdg4GaKz0pdXg24nBns7xF40pe08OsbRq5v1sbGug2u/h7gJWBE0HJ3Es8cEqImxYXs/WGWr3HB+YOVTS4eIVDew7NZyyE5uq8oGfPsVwMMIP33Uxr9+wlH9++VpueP5Kbn38KI8naRvpnMNFy9M3KaywgxFmikTLFicBM591+5ur/YhYD9nJdA0Haar2zSi4wRJgXUPBrCsQH+0dpbWuHO80xeRaa8s5ORQkMo3m1z0css8h9Rgdn9bkIpg/fOQwAB990VlZjbsQJEYrDY6FqSkvrNYA6QmHg8BmEfnEfA9lhXHh0DuNcDg1GJgxisPBsetOp7of6RkhEI6lpTn4PC6ef1Yz9zx3qqQS4na097OswdIUElm3uIajvaPx7M7JxFQ5kUaSlYNTCC9VCeY/PdPJ44d6+dTL1k6oxf+P162mpbaML9yxe9rfbsvhXs5orszI0be0oQKXwKHumav3ZkPcrDODyW02eN0uGiv90wqH3pFQ2vH0zow/20TFo73JE9la68qJxpRT05h4u4eDNKbx+zih0nsSnNLRmPK7p49z9VkL8qqdzZbEPIfBQHhCqG+hSEc4HALuwWrwM69DWWHc6dM3MvXh1jkYmNEJ5rDUtutO11nOib1emoZwALju7AV0Dwd5ugClG9Lluc6hKVEmMB5fnqxk+WBICUVjLE7DbwOwfkktPo8rqWkpHI3x5T8/x1kLq3njxqUT1pX73HzsxWexs2OAO3d1TvneE4d6ueyMiQ2bZsLvcdNaVz6j0z1bnEi5pur8RqYsrJleOPSNhmhI0/HpTJSyrVR7rG8sfp9MJlUeRc9wiKY0BHp1mZfFdeUTch0e3t9N11CQV9u930uVKr+bcFQJRWIMjkVKUzjoadJD2qG+wrro+qbTHIaCLEjDUQeW5K+v8E4bNePYq1sq0xMOV5+1ALdLuGd38forJxIIRzncPTKtcHBq9ySra9M7Zs3gW9IM1fR73Jy/tC6p5vDnZzo50jPKP71w9bR5Aa84bzFtjRV86/4DE4IDth/rZzQ00WeSLiuaKvPmc+geCiICDRX5Fg5lnByc+lDvGwnH74GZcJKyuoYyd86PhaJ0DQXjGvZkHLPjdGHMPSNBGtNserO2tYZdHeONo377VAfVZZ4p3f5KjQrb+TwSjFiaQymalUTkPhG5d/KrEIMrBo6JYbJwiMWs4leTzSipWFRbPq3K3d43SlOVD78nPQdTXYWP85fW8eC+0uhNvf/UMDGFs6ZptNJaW0alz82BJLX+ewLWAzoTlf7itgaeOT44bTLgzQ8doq2xguvOnj7ixO0S3v38lTzdPsBjB8cFzMP7uxGBy1ZmJxwOd49k3EwmHbqGQzRU+PLe1GVhjZ9TQ1Ovzf7REHVpCgfnXujJor6SY25dmsSstMDWSk5NEmCqamkOad6H5y2p5WD3CANjYUZDEf68q5Pr17eURA2lVDhVE0ZCEcvnUIqaA/AR4KP269PAdmBrHsdUVByz0uQoo4GxMNGYZmSfbqktiyfsJNLeN8biJOp0Mq44s4kdHQN5C6HMhOeSxKeD5ZRe3liZ1GHbY2sO6TqkATauaCAa0ynF2J482sf2Y/288/IVKSM5XnvhEpqqfHzr/vG6/w/v72b94lpqs4gdX95YmTISbTZ0DwczmoBkS0Olj77R8ARfTDASZSQUpaEyvd/EuRe6s9AcjtmBGslKYNSUefC5XVNMVoNjESIxTcvnAHDe0jrA0hTvsjsPvur80jYpgWV5ABgJRhkMROKVWgtJOmalbQmvh1X1w8Cm/A+tOPg9bip87inRMfH6Jhn0cF1UW5ZEcxjLOI79ylVNqJL3chLpsKdzEJ/HRVsSk4Azs56OnkCMSp87IzX5wmV1uGRqMtzNDx2iuszD6y5akvL7ZV4377x8Bffv7WL3iUFODgZ48mg/z1+dXaG1VtvvNJ3gny09drRQvqmv8BGN6YRuY46wS1dz8NnhpFlpDr2pNQcRoanKN6WUjXOsdOsMbVjegN/j4r7nTvGrJztYXFeeUV5LsajwW5rNYCDMcLBEfQ4i0pDwahKRFwGlU4QkD9RX+OL1TBx648Wv0p/VtdSU0TMSmtATOhazql5mWo75vKVWgb4H9xU/3+G5ziFWLahKavpoa6rgWN/YtAlovQGlta48rVBJh+oyL+e01k5IhmvvG+XPz3TypouXxWdZqXjLJcup9Ln50p+e4/8ePkQ0prz2wtRCJRmL4sIh92U9ujMwmcwGZ5KTOAlyTKnp+hyc/fQMZ645nBoK4nVLSt9KU7Wf7kn77rezhdPV+Mp9bq44s4kfPHKYB/Z28caNS4veIzodHLOSM7ksSZ8DsA3LjLQNeBT4J+Bv8jmoYlNb7p0wowLipX8zMSs5D5FEu2n3cJBQJJaxcPC6XVx2RiMP7S++38EqPpi8hMjyhkqiMZ1Wa+oZ06xCCC87o5FtR/oYsGe3P3j4MAK843ltaX2/tsLLR190Fvfv7eLb9x/k1RcsTnkOqXCc6fnQHAplVnIEQKJwcDSH+gxMbU1VyQtVpqJ7OEhjpT/lg7qpyk/3pFBW5/9fl4GZ5UPXrcbrFtoaK3jH5W0Zj7UYVNhJiM4EpBiaQ0pxJCIu4C2q+nCBxlMSVJd5GAxMdH5mY1Yaf4iMJeQ9OI1cKmDGTtwTuXJVE3c9e5IjPSMzFv/LhG1H+vjUb5+h2u/h319/XtIIErAcgscHArzonOnag1u0JkSaTDYb9AWVS7JolH79+ha+88BB/rDjOC9Zt4jbthzjZetbMhI0b39eG/WVPo73B3j785ZnPAaH5mo/bpfkvBHRaCjCaChaGM3B1oAThcOgPSvPxL7dVOVLGraciq6hmc1nTVU+dh0fmLDMqTM0U3Z0IucuqeWhG6+hviK9BNZSwNEcTsQ1hxIzK6lqDPhmgcZSMlSXeeM3ioNzE2WicjuaQ2LbUackQDaF1Zywy1yalnqGg7z7R1vpGQ6y+8Qgf//TJ4mmSLbrsev9p3ooJ4tRV1WGQ0pDBgLW4dzFtaxfUsvX79nH+259kkA4yj9cc2ZG+xARXnn+Yt676Yx4qGA2uF3Cwmp/zjUHx7GbyQQkW+ptp3OicBiyJ0TVZen/Nk1VU00/6ZCO+aypyk/PcGhCVFg2wgGs0N25Ihhg3CHtNFOqyeB/kivS+bXuEZHXSiZG4jlOTbknfqM4DIyFqfS5M7rAFk3juHTitjOJ1nFY0VTJ4rpyHsqhcPjug4foHw3x47+5hM+96hx2tA9w73PJ8ynS6VLmrJssHIaCEaKaXQy/iPAvLz+H4WCExw/18vGXrOHMGdp65pNFtWU59zk45pl8Zkc7OJpDzwThYD14qzMwYTRW+hkYm74WWSrSMZ/VVXiJxHRCAynH9JWpcJhrOEX2TgyWqOZg8x7gF0DwdKjKCpZ9b7LPYWAsnPE/qMrvodrvmWB+aO8bo6HSl5YTdTIiwhVnNvHIge6clNIIRWL89ImjvHjdIs5aVM3161tprvbz8xTd55wHfqqKmGVeq8f25AQmx8mfbaP0i5bXc/9Hr2bzRzbxt1euzGofuaIlSQ7LbIhnRxdAOJT73PjcrgnXeTaag6PlTJc0mgxVTU84lFv7HkjQ4gfGwlT5PXnPAyk2ZV4XLoET9j1UksLBrsLqUlXf6VCVFaybYzgYmfAAHhwLZzVbmTzDbO/LPFIpkUtWNjAYiOSkp8D9e7sYGAvz+g1W2Qmv28XLzm3hgb1dEyKsEumw1dyZNJ/pMnDHm5Zkf6E3V/uzdiTnEieHJZeJcIUqneFQU+6JN7EHS7Mr87qmLYSXDOcBn6q17mQGxsKEo0rTDOYz54GYmE8ykOV9ONcQESp9nnhtqZI0K4nIPeksm0/UlHlRJV4VEWzNIYuIgUW1ZXQmPCTb+0azMik5xBvgzFAWOx3u2X2S6jLPhBISm85qJhiJ8ViSSqbH+8co97pnbDzSWDW1o142oZKlyqLaMsbC0QkP19kS9zkUqJH8ZA15KBDOyKQExB/wPSmq5k7GuS7S6TYHTPD/DYyFijKLLgaOdaEYvRwghXAQkTIRaQCaRKQ+IdehDSj9FMNZ4KjViRFL2WYpttSWxZsEqSods9QcltSXs6Daz7Y0G+Ck4uED3Vy6snHCTPGSFY14XMITSQrdnRgYo7WubMY8haZp4t977WKGhW53mA9a4j2Uc+d36B4OUlvuLZjjtLrcO8G3NhiIZGRSAuKZypNDTlPh1GKaybfiaAj9CcLBKkJX+AdlMai0E+Gq/J6i5Gakugrfg5XbsMb+67x+xzyPYHKEQGLZ6ezNSuWcGgoSjsboHg4RjMTinbiyQUTY2NYwa83hWO8ox3rHuHxSVdJyn5u1rTVTSlU4dKRZbrvZjn9PNLs4pT+y9TmUEgtrrAfbdCWls8Wywxfut6kp80yYlTtNZTKhIUWhymSMm8/S0xwSzUrDwcwF2FzFCWcthA9qOpIKB1X9L1VdAXxEVVeq6gr7dZ6qzmvhENccEkwG2TbcaKktQ9WyyTphrLPRHMByzHb0j80qWuZhu7Pc5dNUJb1gaR072gemDWk93j9GaxoVVRurfIQiMYYSIk16R0K4BaqzcMaXGvG+H1mUjkhGoRLgHGrKp5qVMp2VV5d5cMnUWmSpSNfx7kzGEh3SI6FIVsEcc5Fmu+PkgpoSEw4JxESkzvlgm5jel78hFZ/KhIqIYDUIGQpGsnZIgxXO2p6YADcLHL/D1lloD48d7KG52s+ZC6b27z27pcZqr9g3UfgEI1aZ5XQ0h6ZpzA19oyEqvZJR6YxSJR4KmkWMfzK6h0MzzqZzSU3ZxJDtoSzMSi6XUF/hm7Y5VjK6hoK4XTJjlnO514qo6h8b3/dwIBKfUc93nJ4n6bQlzgfpCId3q2q/80FV+4B3521EJYCTuj4a7+FqZ45m4ZAe79cbGM9xmKXmcHZLNRU+N1tn4XfY0THA+Uvrpn1Qr7LbK+49OTEiygndbE2jUc94OefxG7t3JESBAnHyTk25B49LMnLEzkT3cLAgOQ4ONZOSPUeDkawcn3UV3oyqBVulM3wz2tFFZEpE1XDw9BEOTqtYx/dQaNIRDu7EBDgRcWN1hZu3ODfIqK05OEk4VVnYOltqxktotPeNUlfhnfXF7XG7OH9pHU8e7c/q+8PBCIe6R1jXOn39xFV2e8XJ4bLH0wxjhSSaw0iYKu/c1xrAenA1VProzZHmEAhHGQpE0q42mgtqyr0EI7F4D+iRUDQ+McoEq1Bl+malgbH0GwoltsuMRGMEI7HTxqzkTEYj0eK0B05HOPwZ+JmIXCsi1wI/tZfNW+KaQ8i+aWwNIptZVU25h3Kvm07brDRbf4PDuYtr2dM5lHFmKsDuE4OoWv2ep6OmzMuimjL2nZzYsCed7GgHx7GaGM7aOxqi2jc/hANYfodcaQ7OfgppVnJMSI5paSwUpSKLB299pS8jh/TgWPrmq0rfuHCI34eniXBwHPLFCuBIRzjcCNwHvNd+3QN8LJ+DKjYVcc3BmVFZF2c26p2IWAlTgwGO9o6ypG52/gaHdYtrCUVj7MsiGe4Zu23iusXJK6+vWlg1Zd+OcEinj3ZDpQ8Rq7OZQ/9oiKp5JBwaq3w5c0g7GlYhHdJxDTkYJRyNEYrGqMiiQ1p9hTcj4TAUTL/aQJWdkOp8D6z+yqcDLzpnEZ++fi0fum5VUY6fToZ0DPgBcJOqvk5Vv62q06fPzhPKvC5Exs1Kzswl2xnLotoyjvWOcqRndFoHcDY4D/ZdHZlXMnmmY5CmKn/KftirFlRb7UATIpaOD4zRVOVLq8Wix+2ivsIXL3Ueiyl9o2Gq54lZCaCh0j+lKVS2jEfwFG6W6Ex2RuxqsED2msNIOO1s8Uwc31V+T3xydrppDm6X8DdXrJhVkcjZkE6G9CuwWoP+2f58voj8Ps/jKioiQoXXHb8Ynb/Z2GPBEg5OaGiuhMPyhgqq/R52dgzMvPEkdh0f4JzWmpRRQ6sXVhEIx+K9fiH9HAeHpipf/KE3FIgQjen80hwqs2t0Mx2FrKvkUJ7gW3MmQtn6HELRWFzAzMTgWDh9s5LfE7//4r6/00Q4FJt0zEr/AlwM9AOo6nZgRf6GVBpU+D2Mha2L0blxsr0oz00w3+RKOLhcwtrWGp45nplwiMaUg10jrJmm/3MiZ9jjPJjQ7vNEmjkODo2V4+WcnVDH+eRzaKz0MRSMxB26s8H5nWYqKZFLKm1BMBKMjmsOWQiHTBLhVJWhQPptL6v87rhQGDHCoaCkIxzCqjr5CVQc93kBqfQlag7OrCq7izIx0cyJBMoF6xbXsvvEYMr+C5Np7xslFI2xsjl18boVdnG7Q12WcFBVKwEuE82hery+kmN+qZpHZXGcvhSZROoko2soSJXfk5bJLlck+tZG49pxdqGskN7vEAjHiMQ07UzsRIf08CzNu4bMSOdX3iUib8YKaV0FfAB4JL/DKj7lPk+CQ9qxdWZ3465aUMVrLljMy89vxe/J3c1/1sJqAuEY7X2jaXeGO9BlRSCd0ZxaSDVW+qgu83DI1hwGxyKMhKJp5Tgk7sMp0+38nW9mJbCa3s/kpB8KhPnKn/cwHIzwiZesYcGkbniFLp0B49dzolmpMhuzUmX6msNgvGdE+mal0VCUWEwZs+/D8gIK0NOZdP5D/wDcBASxwljvBD6fz0GVApU+9wSHtEj2F6WI8B9vPD+Ho7M408lHODmctnA4aGsCMwkHEWFlU2VcOHRkEMbqUF/hYzAQIRyNxR8c880hDellSd/0m2f4/dPHAUt7+/l7Lpvg8yl06QwY1xJGQuNmpfIsfQ6QnnCIJ5SmG62UUK0gYJvvCqldnc6kE600qqo3AdcCV6vqTaqa+87qJUaF3xPXGEaCUSp9npIr+7DK9gtkEs56oGuYhkpfWrHTKxKEQyY5Dg5O34a+0VD8wTGvNIcqp75S6odi+1CM3z99nA9ccyZffPU6thzu49EDE0uip9M2M9fENYdgJCFcO4topbhZKR3NIbOGQvFSNsEogbCV01Pmnd+NfkqFdKKVNorITmAHsFNEnhaRi/I/tOJS4XUzmuAIK1YKeyqqy7y01k5NVkvFgVMjrEyzWU5bUyXHB8YIhKPxIn+ZmJXi5oaRML0jYXxuF2Wl9zNmzbhZKfVD8f72MD6Pi3devoLXXriE2nIvv9zWPmGb7uFgQXpHJ1LmcUJZEzSHLGblteVeRKAvjeJ7TrmOdAv8OVUJhhMc/0ZzKAzpiOCbgfepapuqtgHvB76f11GVAOU+N0E7+3gklF3NmUKwamH1lBpIqTjQNTyjSclhRVMlqnCkZ5SO/gA+t4umDBrROFEsvSMh+kZC1Fd6S077mg2OUzWxPtFkVJVtJ6NctbqZ+korR+TaNQu4b88pIlHr+gpHY/SPhguuObhcQoXPmgSNBrMPZfW4XdSUpZcI52Rjpxut5CTljYXGNQd/gfpdnO6k8ytHVfVB54OqPgTkrv1VieL3uOKtMsdC0ZKdraxaUMX+U8NpRSwNjIXpGQnNGKnksLLJEiKHuoc53j9GS11ZRk1HGhLMLr2joXnRAS4Rt0uo9num9BtP5LnOIXoDygvWLowv27RmAX2jYZ49YSUw9hahdIZDhc8yn47ZD95sI/KsLOk0NIe4Qzo94eDcd4FIlGA4it/jmlcTjFImHeFwv4h8W0Q2ichVIvI/wGYRuVBELpzNwUXkn0RERaQpYdkmEdkuIrtE5P7Z7H82+D2uuOYQjMRK1s65emE1wUgs3isiFcd6rW2WNaRXwqOtydruYPcI7X2jGeU4QILmMGppDvOhA9xkasq9E/oNTMbpqPe8hKZKFy2vB+Apu3Ci03+5ucBmJbD8DqOhSHwilO2svL7Sl1Zl1rjmkGZvFOe+C4SjBMKlO0mbj6TzHzrP/vsvk5ZfgJXvcE02BxaRpcALgaMJy+qA/wFerKpHRWRBNvvOBX6vO27jLOWLss32HxzpmTmc1REOS9MUDtVlXpqr/RzqGuFA1wjXr2/JaGx1ThTLiOWQXtNSA8yvWIaacm/KPtJPHOqloUwm9PBorS1jYY2fJ4/28fbntXFy0PpNFtYUvm5/uV0JIBCJ4vO4sm5HWV/hi59HKoYCYdwuSdu3EdccwjEC4dKdpM1HZhQOqnp1no79n1gF/H6XsOzNwK9V9ah97FN5OvaMlNmag6oSiERLtjVhW6P10DnSMwI0p9zWKYWxNINmQyuaKtl2pI+BsXDavgoHn8dFtd9j+RxGw/GolvlEbbknpc9hR0c/K2snPtBEhAuW1sc1hxN2n4x0ChrmGsu3FiUYjlE2C1t+fYWPPZ0z+76ciqzpmoYmaA6RaE7zhAypSfrEE5GXAztU9Yj9+Z+B1wJHgA+q6qFsDyoirwQ6VPXpSRfJasArIpuBauC/VPVHSfZxA3ADwMKFC9m8eXO2w5mWjmMhVOHu+zbTOzBGWWQk58cYHh6e9T5VFZ8bHtq+h6XBwym3fezZIBUeeOqJh9Pef1k4yMFuO2rrxAE2bz6S0fjK3VGePXiMvpEog10nGPaEcv47FpPQcIBTo7Fpz2kkrBzrHeP85TplfU0kxNHeML+/8z4ePxpGgGe3PcaeAjeSDwyPMTIEMtaPS6NZ/2+Ge4N0DUXi3092be8/EsCr0/9e09E1apl2t+/cRXtXlGgo/e8Wmlzcz6VEqunwF4FLAUTkeuAtwJuwzEnfAl6UascicjewaJpVNwGfxDIpTTeei7ByKsqBR0XkMVXdO3lDVf0O8B2ADRs26KZNm1INJ2P2uw/Cvt1c8rwrcG99kKWtDWzadH5Oj7F582ZyMe6VTz9AtKKcTZs2ptzuB4eeYMWCIJs2XZn2vntr2nng508D8MYXX5Fxy8LWXQ8zGouhDHLB2lVUhY/k5JxLhdu7nqZzf/e05/TYwR645zFWNZVNWV+xvJef73mU8qVrKRvsZEFNF9deky8lPTk/PPQE3cMhGpoqqRntz/p/s0v3c+fhPVx6+ZWUed1s3ryZK658Ph73RG3kx4e3sEACaV+DXUNBeOBu2s5YxYHwKcLeEJs2XZHVGPNNru7nUiGVHqmq6ng5XwPcrKrbVPV7zGS/sL58naqum/wCDmIV7ntaRA4DS4AnRWQR0A7cqaojqtoNPMC4z6OgOI65oG3rLOXwueWNFRzuSc8hnYlJCeDF6yz5vums5qx62TZU+jhwaiT+fr5Rm8Ihvc8OMV5SPfXaWbe4BrdL2NHeT+dggEUZOvtzRbnPzVg4Omt7vhOJ1m9HLN3fHmbdZ+7kz890TthuMJB+RVZINCvFLN+fMSsVjFRXg4hIlYi4sGby9ySsy9o4qqo7VXVBQt5EO3ChqnZi+R+uEBGPiFQAlwC7sz3WbPDbjrCgHUJXqg5pgLbGSo72jk7ovTAZVaW9b4ylDZk9hCp8Hp745LV86y3Z5T3WV/gYsyNh5lsoK1jx+iOhaDxnIZEDXSNU+tzU+6eaiip8HlYvrGb7sX4O94zkrENgppR53ATCUcbC0VnVLHL8Sb0jISLRGD/fEyIQjvGlP028fTOpyAqJDmlLgPmNQ7pgpPqlv4bVx2ErsFtVtwKIyAXAiXwMRlV3Y/WN2AE8AXxPVZ/Jx7FmwtEUAuEYgUhpX5TLGisIRWJ0pogW6RoKEozE0o5USmRBTVnWwjEx67cY0Tj5xsngdSr4JnKwe4QVzZVJna/nL63jwX3dHOsd45zW6Vu25psynzseJuqfjXCodDSHEE8c7mUkDJeubOBwz2i89Ao4vRzSFw5etwu3SwhETChroUn6xFPV/wOuAv4GeGnCqk7gnbkagK1BdCd8/qqqrrXNUF/L1XEyxYmKCISjhCKxko6SaLNDWA8n9F6YzLE+6wbN1Kw0WxoTTEnLGwt77EJQbdf+cVpYJnKwazieSDgdV60et86ua03esjWfWJqDNQGazYO3PiGnZcuhPgT44LWrAXjyaF98u0y6wI2P0UUgHLPzjUr3PpxvpJwOq2qHqj5ltwp1lp1wQk3nM46m4GR0lnJ8tZPU1t43lnSbTjtcsiWD2ki5ILHR0Xy8sRMLwyUSjsY43j+WUiA+f3UTVX4PHpewfklxhEO5z6oEEAxHZxfKGi+yGOaZ4wMsqhQ2tNXjc7vi3QpjMWU4FEm7IqtDmTdBuylh3998ozSD90sAx/HlxLCXsiNsYU0ZIlaP52Q4JqeFWTiVZ8P5y+oKerxC4xRkHJ6kOXQOBIgpli8hiUJX4fNw70euotLnKVoDmzKPm0hMGQ5GZiW868pts9JIiF0dAyyvceF1u1jSUM5RO1hiKBhBNf2ie/Exet0JmoMRDoXCCIckOJqDE4lSyj4Hn8dFc5V/gm13MqcGA/g8rnjXrkJR4fPw9suWs7ZINvV8Ux2vGjpRc3D6XyyuqyCS3NqXVQRYLnEEwsBoeFYPXp/HRZXfw/6uYY4PBLhykSUsljdUcMQRDk4vhwx8DtYYXeM+hxKepM03UgoHEXEDu1R1TYHGUzI46uvAHNAcwOqzcLw/uUO6czDAopqyohQt++wr1xX8mIVi3Kw0sYRGh23iW1xfzpGOgg8rbcrsKqxDs9QcwDItPbTPch+21Vj3z/LGSrYc7ov3job0eznEx+h1EzS1lQrOTD6HKLBHRJYVaDwlQ3xG5QiHEr8oF9eVp9QcOgcCLKwpfNXP+Y7TqWw4MFE4OP+LliKUxMiERD/DbNtv1lf44r0tltnCYVlDBcPBCH2j4biJNpNoJbDuvaFAhJiWtu9vvpHOL12P1Uf6HhH5vfPK98CKjaM5OEXVSt0R1lpXRkf/GKrT5zqcHAzMy1DSYhMXDpM0h5NDARrs/g2lTGJb0NmEsgKcY0dcrWiqpNJuB+vUizo1FMi4IqtDmdc1ZyZp84l0/kufzvsoShAndHWuXJStdeUEIzF6R0I0Tmoao6qcHAxy3dlGOOSayiTC4dRgkOYCN+/JhkRz6Wxn5VetbuKnTxzl4rYGwCpV3mz3qOgaCsbDfTPWHDzueOb1bAWYIX3S6SF9P3AY8NrvtwBP5nlcRWcuOaRhvLfzdH6HwUCEsXDUaA55wOt24fe4pvgcuoaDLJgDZrzESc9sux1ed/ZCPvWys/nky86OL3MEZNdQMK6FZ+NzcLrMzSbc1pAZ6fSQfjfwS+Db9qLFwG/zOKaSwGcXDBu1G6/73KV9US62hUPHNH6HeL+AErd/z1Wq/J45qzmU+8av69mG03rcLv72ypXUJuQxTNAc4l3gMjuO3zveeMtoDoUjnSfe+4HLgUEAVd0HFK0JT6Hw2sLASW7ylrhwGNcckguHRUZzyAsVfvcEzUFV6RoO0jwHNIcq//iDvDKL/tEzUen3UOFzW5pDIILf48q42kCidmM0h8KRjggPqmrICYEUEQ9WB7h5jdsliCRoDp7S7ltbX+HF73FxYppEOCc72kQr5YcKrydeXBAsM14oEpsTmkPiLD5fiXjN1X66hoNU+NwZ+xtgsl/EaA6FIt0e0p8EykXkBcAvgD/kd1ilgdftYiQ0NzQHEWFRbRmdg8Ep64rZhvJ0oMznZiw8XpW1Z9j6HzQWoSd0phRCODRV+eOaQ6aRSjDRUW6EQ+FI54n3caAL2Am8B/gj8Kl8DqpU8LldjM0R4QDWw//kwFSH9MnBIHUVXnNj5YkKr5tAaFxz6LeDGOrmQInyRCe0Uwok1zQ7wiHDiqwOE8xKJR4YMp9IR4xfDdyiqt/N92BKDa9b4jf6XBAOi2rK2H6sf8ryzsFAwWsqnU6U+9ycGhoXyv12ZM1c6F/hSmhLOttopWQ0V/t57FAPlX5PxnWVwGgOxSKdJ97bsLq2PSYiXxWRl4tIfb4HVgp43C6cnLJSj1YCbLNSYEoi3MnBgIlUyiPlXndcw4Txbmh1GVYfLTZVefQ59I+G6RkJZlxXCSZmbpd6GZv5RDp5Dm9X1dVYrUKPAf+NZWaa9yQKBG+JO6TBMiuFIrH4w8nh5GCARcYZnTfKfVbVUIc++/efC5pDIhX5MivZ4azHescyDmOFieGrpZ5vNJ+Y8T8lIm8BrgTOBbqBbwIP5nlcJYHXLQnvS/+idEJVOwcD8c5ckWiMrqGgcUbnkXKve0K00sBoCJdkHs9fbPLV0CoxaivTXg4wOZTVaA6FIp2r92vAAeBbwH2qejifAyolEgWCx1X6msOiWusm7BwMcHaLVSK7ezhETE2kUj4p97njIc9gaQ615d4J9vzTmcRM8dpshENCboPRHApHOmalJuBdQBnwRRF5QkR+nPeRlQCOcPC5XUUpdZ0pjgBIjFjqNAlweafcbkYTi1m+nr7R0JyIVHJw51mItTVVxt8vy6KHeaLmUOoFMOcT6ZiVaoBlwHKgDagFYqm+M19wzEqJ5qVSxmkc4wgEMDkOhcCpbBqMxCj3ue14/rnjjH7k49dM8VPlkkQn9BnNyXtqJyNROMyFSdp8IR2z0kMJr2+qant+h1Q6OJqDd47MVnweF01VvrhAgMS6SsYhnS+caJrRUIRyn5vhQJjqIrX9zIaFNWUFmzysSNAi0sXkNhSHGa9gVV0PICKZi/w5Tlw4zAFntMPCmrJ4uQywSme4XUJTpREO+cLRHByn9EgwGo/QMVhcv76FO3aemNA/Il0cR3nTHChHMp9IpyrrOhF5CtgFPCsi20Rk/vZ9TMDRGOZCjoPDopqJJTRODgZZUO03ztE84mgOAVs4DAcjEwraGeAbb7qAA198aVbfrbAFysvPa8nlkAwzkI7u+x3gw6p6H4CIbLKXPS9/wyoNfHPM5wBWWe6nErKkTQe4/OM8vEZDicLBhFwmImIVssyGpQ0V3PGBK1izqCa3gzKkJJ0pcaUjGABUdTOQueFwDuJxzT2z0qKaMnpHQgQj1oOqczBgIpXyjGP2CEZiqCrDwUjeitidrpzTWpv3qCrDRNJ56h0UkU+LSJv9+hRwMN8DKwUcs9JcEw5gNZsBR3Mwttp84rOvk1AkRjASIxpTquZYApzBMJl0nnrvApqBXwO/Apy8h3lPPJR1jkQrwXi3t87BAKOhCEOBiKmrlGec2PtgJBrvCJevOkUGQ6FIegWLSBnwd8CZWOW6/0lV8xcMXYL44klwc0edjZfQGAjQaJfQMGal/JKoOQwHLOGQrwqnBkOhSHUF/xAIY9VReglwNvChAoypZPDEHdJzR3NwBMHJQSMcCoUvrjnExjUHY1YyzHFSXcFrVfVcABG5GXiiMEMqHeZinkNNuYcyr4vOgUA8S3dxfXmRRzW/cTTMCcLBmJUMc5xUV3DchKSqkdMxbd03B4WDiNi5DgEq/R5EoKXWCId84hSDC0VijNjCwUQrGeY6qa7g80Rk0H4vWD2kB+33qqrzPug4XnhvDvRySGRhTRknBwOUe90sqPbHzR6G/OB3W6GsoUgsniVdbjqWGeY4SYWDqp72V7cTV+2aY1rTotoyth3pw+t2sbjOaA35JtHnELSb/ph6QIa5jrmCU+Dc9OsW1xZ5JJmxsqmKjv4xdnYM0NZ4WuQrFpXEaKWAnXxoeh0b5jrGMJqCv75kGWtbath0VnOxh5IR57TWoApDgQjnLplbgm0u4nYJHpcQikbj7UJNxzLDXMcIhxTUVfi4es2CYg8jY85ZPO4OWm+EQ0HweVwEw7F48T3Tscww1zHCYR6yqKaMt166nKO9o3POJDZX8XlchKIxguEoIqZjmWHuY4TDPERE+PyrTouq6iWD3+OyfQ4x/J650VbWYEhFUaY3IvIZEekQke3266UJ69aLyKMisktEdtplPAyGksbncdnRSlHjjDbMC4qpOfynqv574gIR8QC3AG9V1adFpJGEZDyDoVTxuW3NIRwzJiXDvKDUzEovBHao6tMAqtpT5PEYDGnh87gJRmJ43EZzMMwPRFULf1CRzwDvAAaBrVgVX/tE5EPARcACrDLht6nqV5Ls4wbgBoCFCxdedNttt+V/4DlmeHiYqqrTqzX3fD3nzz06RoVH8Hvg5EiML1xREV83X885Feac5wZXX331NlXdMN26vAkHEbkbWDTNqpuAx4BuQIHPAy2q+i4R+QjwfmAjMArcA3xKVe9JdawNGzbo1q1bczn8grB582Y2bdpU7GEUlPl6zm/49qMIVvJb/2iI3/39FfF18/WcU2HOeW4gIkmFQ97MSqp6XTrbich3gdvtj+3AA6raba/7I3AhlpAwGEoWv8cVr8jqN2YlwzygWNFKLQkfXw08Y7+/EzhXRCps5/RVwLOFHp/BkCk+t5UEF4zEjM/BMC8olkP6KyJyPpZZ6TDwHgDb7/AfwBZ73R9V9Y4ijdFgSBuv20UkFiMWVhZUm57dhrlPUYSDqr41xbpbsMJZDYY5g8ctRKKKokZzMMwLSi2U1WCYk3jdLsKxGJGomnLdhnmBuYoNhhzgdgnRqBKOxkxzJcO8wFzFBkMO8LqFcEwJRWJ4XOa2Msx9zFVsMOQAj8tFJBojElO8blN0zzD3McLBYMgBHrcQiVlmJaf3uMEwlzFXscGQA7xuF5GoEo4qHiMcDPMAcxUbDDnA7ZJ4/2ifMSsZ5gFGOBgMOcDrEpwyZUZzMMwHzFVsMOSARIHgcRnNwTD3McLBYMgBngRTkslzMMwHzFVsMOSARG3B5DkY5gPmKjYYckCiQDB5Dob5gBEOBkMOSBQIJs/BMB8wV7HBkAMmOKSN5mCYBxjhYDDkALfLaA6G+YW5ig2GHDDRrGQ0B8PcxwgHgyEHTHRIm9vKMPcxV7HBkAMStQUTymqYD5ir2GDIASaU1TDfMMLBYMgBbhPKaphnmKvYYMgBXpcJZTXML4xwMBhywITaSkZzMMwDzFVsMOSACQ5pIxwM8wBzFRsMOcDtMiW7DfMLIxwMhhyQKBBMyW7DfMBcxQZDDvCaZj+GeYYRDgZDDkh0SHuN5mCYB5ir2GDIAYnagtdkSBvmAeYqNhhygCnZbZhvGOFgMOQA74Q2oUY4GOY+nmIPwGCYDzRV+fnbK1bwho1LETHCwTD3McLBYMgBLpfwqevXFnsYBkPOMGYlg8FgMEzBCAeDwWAwTMEIB4PBYDBMwQgHg8FgMEzBCAeDwWAwTMEIB4PBYDBMwQgHg8FgMEyhaMJBRP5BRJ4TkV0i8pWE5etF5FF7+U4RKSvWGA0Gg+F0pShJcCJyNfBK4DxVDYrIAnu5B7gFeKuqPi0ijUC4GGM0GAyG05liZUi/F/iSqgYBVPWUvfyFwA5Vfdpe3lOk8RkMBsNpTbGEw2rgShH5IhAAPqKqW+zlKiJ3As3Abar6lel2ICI3ADfYH4dFZE8Bxp1rmoDuYg+iwJhzPj0w5zw3WJ5sRd6Eg4jcDSyaZtVN9nEbgEuBjcDPRWSlvfwKe9kocI+IbFPVeybvRFW/A3wnT8MvCCKyVVU3FHschcSc8+mBOee5T96Eg6pel2ydiLwX+LWqKvCEiMSwpG478ICqdtvb/RG4EJgiHAwGg8GQP4oVrfRb4GoAEVkN+LDUsTuBc0WkwnZOXwU8W6QxGgwGw2lLsXwO/wf8n4g8A4SAt9taRJ+I/AewBVDgj6p6R5HGWAjmtFksS8w5nx6Yc57jiPVMNhgMBoNhHJMhbTAYDIYpGOFgMBgMhikY4VBgROSfRERFpMn+LCLydRHZLyI7ROTChG3fLiL77Nfbizfq7BCRr9olUnaIyG9EpC5h3Sfsc94jIi9KWP5ie9l+Efl4UQaeQ+bb+QCIyFIRuU9EnrXL3HzQXt4gInfZ1+tdIlJvL096jc81RMQtIk+JyO325xUi8rh9bj8TEZ+93G9/3m+vbyvqwLNBVc2rQC9gKVZE1hGgyV72UuBPgGDlfTxuL28ADtp/6+339cU+hwzP94WAx37/ZeDL9vu1wNOAH1gBHADc9usAsBIrgu1pYG2xz2MW5z+vzifhvFqAC+331cBe+3/6FeDj9vKPJ/y/p73G5+IL+DDwE+B2+/PPgb+y338LeK/9/n3At+z3fwX8rNhjz/RlNIfC8p/Ax7AisRxeCfxILR4D6kSkBXgRcJeq9qpqH3AX8OKCj3gWqOpfVDVif3wMWGK/fyVW9ntQVQ8B+4GL7dd+VT2oqiHgNnvbucp8Ox8AVPWEqj5pvx8CdgOLsc7th/ZmPwReZb9Pdo3PKURkCfAy4Hv2ZwGuAX5pbzL5nJ3f4pfAtfb2cwYjHAqEiLwS6FC7blQCi4FjCZ/b7WXJls9V3oU1e4TT55zn2/lMwTaXXAA8DixU1RP2qk5gof1+vvwOX8Oa3MXsz41Af8IEKPG84udsrx+wt58zFCvPYV4yQ8mQT2KZWeYVqc5ZVX9nb3MTEAFuLeTYDPlFRKqAXwEfUtXBxImxqqqIzJs4eRG5HjilqttEZFORh1MQjHDIIZqkZIiInItlW3/avoGWAE+KyMVAB5YvwmGJvawD2DRp+eacD3qWJDtnBxF5B3A9cK3aBliSnzMpls9FUp3nnEZEvFiC4VZV/bW9+KSItKjqCdts5FRbng+/w+XAK0TkpUAZUAP8F5aJzGNrB4nn5Zxzu13toRaYW1Wmi+30OB1fwGHGHdIvY6Kz7gl7eQNwCMsZXW+/byj22DM8zxdjlT9pnrT8HCY6pA9iOW899vsVjDtwzyn2eczi/OfV+SSclwA/Ar42aflXmeiQ/or9ftprfK6+sCZtjkP6F0x0SL/Pfv9+Jjqkf17scWf6MppD8fkjVjTHfqxKtO8EUNVeEfk8VikRgM+pam9xhpg138QSAHfZGtNjqvp3qrpLRH6OJTgiwPtVNQogIn+PFdHlBv5PVXcVZ+izR1Uj8+l8ErgceCuwU0S228s+CXwJq8Ly32BF5L3BXjftNT5PuBG4TUS+ADwF3Gwvvxn4sYjsB3qxBMScwpTPMBgMBsMUTLSSwWAwGKZghIPBYDAYpmCEg8FgMBimYISDwWAwGKZghIPBYDAYpmBCWQ2nNSISBXYmLHqVqh4u0nAMhpLBhLIaTmtEZFhVq5KsE6x7JDbdeoNhPmPMSgZDAiLSZvdf+BHwDLBURD4qIlvsXgSfTdj2JhHZKyIPichPReQj9vLNIrLBft8kIoft9267x4Wzr/fYyzfZ3/ml3f/iVqeCp4hsFJFHRORpEXlCRKpF5AEROT9hHA+JyHmF+o0MpwfGrGQ43SlPyPI9BPwjsAp4u6o+JiIvtD9fjFX+4fci8nxgBCvr9Xys++hJYNsMx/obYEBVN4qIH3hYRP5ir7sAq6zIceBh4HIReQL4GfBGVd0iIjXAGFb27TuAD4nIaqBMp1b7NRhmhREOhtOdMVU93/lgl6A+olbfAbAq6b4QqzQCQBWWsKgGfqOqo/b3fp/GsV4IrBeR19mfa+19hbDqDbXb+9oOtGGVeT6hqlsAVHXQXv8L4NMi8lGsUug/yPCcDYYZMcLBYJjKSMJ7Af5NVb+duIGIfCjF9yOMm2zLJu3rH1T1zkn72gQEExZFSXFvquqoiNyF1VDmDcBFKcZiMGSF8TkYDKm5E3iX3bsAEVksIguAB4BXiUi5iFQDL0/4zmHGH9ivm7Sv99rlrhGR1SJSmeLYe4AWEdlob19tl38GqxvZ14EtanUKNBhyitEcDIYUqOpfRORs4FHbRzwMvEVVnxSRn2GV4T7FePVcgH/Hqk56A3BHwvLvYZmLnrQdzl2Mt5Wc7tghEXkj8A0RKcfyN1wHDKvVdGYQ+H5uztRgmIgJZTUYcoCIfAbrof3vBTpeK1bzpzUm1NaQD4xZyWCYY4jI27B6Nt9kBIMhXxjNwWAwGAxTMJqDwWAwGKZghIPBYDAYpmCEg8FgMBimYISDwWAwGKZghIPBYDAYpvD/AQ/A+m3RFZZCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "import math\n",
    "n=np.arange(24)\n",
    "fs=1000\n",
    "x=np.exp(2*np.pi*0.4*n*1j)+np.exp(2*np.pi*0.425*n*1j)+0.2/math.sqrt(2)*np.random.normal(size=n.shape)+np.random.normal(size=n.shape)*1.0j\n",
    "plt.psd(x,NFFT=150,Fs=fs,window=mlab.window_none,noverlap=75, pad_to=512,scale_by_freq=True)\n",
    "plt.title('Fres=41.6888mHz')\n",
    "plt.ylabel('Power Spectrum(dB)')\n",
    "plt.xlabel('Frequency')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "637c56f7",
   "metadata": {},
   "source": [
    "使用子空间方法来解析两个紧密间隔的峰。在此示例中，使用 MUSIC 方法。估计自相关矩阵并将自相关矩阵输入到pmusic中。指定具有两个正弦分量的模型。绘制结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "835fa33c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<function matplotlib.pyplot.show(close=None, block=None)>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAm40lEQVR4nO3de5xcdX3/8dd7d5NsIBduAUwCBAW1iBd0sSItUkF/GrkIYtHWC0WlXqhQUbxgW3uzVizya9GfRmhFpUBBYlFiuSiIWLkk3AQCCIoSAmQJSC4k2ct8fn+cM7snu+fMzmZnzsxu3s/HYx+ZOfOdcz5zZjKf+V7O96uIwMzMrN10tDoAMzOzPE5QZmbWlpygzMysLTlBmZlZW3KCMjOztuQEZWZmbckJykon6SRJN7U6jqlO0h9KeqDVcYyXpM9IOr/VcVjrOUFthyQ9ImmTpA2SnpT0TUmzWh1XO0jPzZFtEEP1/an+nVfH80LSftX7EfHTiHhRk2L8pqR/aMa+I+LzEfH+bXmupBvS8/DyEduXptsPT++Pil/SorRMV3r/DyT9r6RnJT0t6WeSDk4fG/UjS9KfSFqevl+PS/qhpD/YltdhCSeo7dfRETELeCXQA3y2xfFMCtUvrxIcHRGzMn+nlnTcqeBB4D3VO5J2BQ4BeuvdgaQ5wA+AfwN2ARYAfwtsKSj/MeBc4PPAHsDewFeBY7flBVjCCWo7FxGPAT8EDlTiy5LWSFon6ReSDgSQNEPSlyT9Nq11fU3SzPSxvF+TQ7/mJe0q6cp0n7cCLxhR9rWSbkt/qd4m6bWZx06S9CtJ6yX9WtKfZrb/TNJ56fPul3RE5nlzJV2Q/pJ9TNI/SOrMPP4BSSvT/d4n6ZWSvk3yxfL99FfwmZlf1e+T9Fvgx5IOl7RqxGsYqnlJ+pykyyR9J93/LyS9UNKn03P7qKQ3bsv7JWk/ST9JX/NTki5Nt9+YFrkrjf3EkXGmMX5C0t2SNqbnZ4/0l/56SddJ2jlT/jJJT6THulHSS9LtpwB/CpyZHuv76fb5kr4rqTd9rz5a8Bp+P91v9v04TtLdmfP3nbHiqOEi4MTM/t8JLAX66jnHqRcCRMTFETEYEZsi4pqIuDvn9cwF/g74SERcEREbI6I/Ir4fEZ8YxzFtBCeo7ZykvYDFwB3AG4HDSP5zzgX+GFibFv1Cuv0VwH4kvyj/us7DfAXYDDwPODn9qx5/F+Aq4F+BXYFzgKvSpLZjuv3NETEbeC1wZ2a/vw88DOwG/A1wRbo/gG8CA2msB6Wv7f3pMd8OfI7kV/Yc4BhgbUS8G/gtw7WXL2aO9Trg94D/U+drPhr4NrAzybm9muT/2wKSL7Ov17mfkf4euCbd70KSX/hExGHp4y9PY7+04PlvA95A8l4eTfLj5DPAvDS+bFL5IbA/sDtwO8kXPxGxJL39xfRYR0vqAL4P3JW+xiOA0yWNOl8RcQuwEXh9ZvOfAP9ZEHNuHDWsBu4jec8heZ+/NcZzRnoQGJR0oaQ3ZxN3jkOAbpIkaA3kBLX9+p6k3wE3AT8haZroB2YDLwYUESsj4nFJAk4B/jIino6I9Wn5d4x1kPRX7NuAv05/Wd4DXJgp8hbglxHx7YgYiIiLgftJvjwBKiS1u5kR8XhE3Jt57hrg3PTX6qXAA8BbJO1BknRPT4+5BvhyJt73k3y53haJhyLiN2O8lM+l+9o01mtO/TQiro6IAeAykgTwhYjoBy4BFknaqcbzvyfpd5m/D6Tb+4F9gPkRsTkixjvY5N8i4sm05vxT4JaIuCMiNpN8wR5ULRgR/x4R6yNiC0lCf3laW8hzMDAvIv4uIvoi4lfANyj+jFxMUrNB0myS9+vivILjjKPqW8B7JL0Y2Ckifj5G+ZHHXAf8ARDp6+hV0gqwR07xXYGn0vfaGsgJavv11ojYKSL2iYgPp00YPwbOI6nxrJG0RElb/DxgB2BF9QsT+J90+1jmAV3Ao5lt2WQwf8T96uMLImIjcCLwQeBxSVelXzhVj8XWsx3/Jt3fPsC09DnVeL9O8gscYC+Smtd4PDp2ka08mbm9ieQLbDBzH6DWwJTq+1P9+0a6/UxAwK2S7pV0co191BPXyPuzIPlhIekLkh6WtA54JC2zW8F+9wHmZ5MqSc0s7wsdktrS8ZJmAMcDt+f9SNiGOKquIKmhnUpSkx1pgOQzkjWN5AdRBSD9gXZSRCwEDiT5bJ2bs6+1wG4qr39yu+EEZVuJiH+NiFcBB5A0A30CeIrky+slmS/MuekgC0iaa3ao7kPSnpld9pJ8GeyV2bZ35vZqki83Rjz+WBrP1RHxBpLmwftJfs1WLUhrd9nnrSZJJluA3TLxzomIat/Fo4zoB8uegjq2j3y9ndSXrCcsIp6IiA9ExHzgz4GvKjNyr4H+hKSD/0iS5t5F6fbq+R55nh4Ffj0iqc6OiMV5O4+I+0h+ULyZ2s17Y8WRKyKeI2ka/BD5Ceq3mX1V7Qs8GhGVnP3dT9JsfGDOvn5O8nl7a62YbPycoGyIpIPTDuxpJF/Cm4FK+h/2G8CXJe2ell2Q6V+4C3iJpFdI6iZphgEgrTVcAXxO0g6SDgDemznsMuCFSobodkk6kSQ5/iDtwD827YvaAmwg/XWb2h34qKRpab/S7wHLIuJxkn6af5E0R1KHpBdIel36vPOBj0t6lRL7SaomySeB549xqh4EuiW9JT1XnwVmjHmCG0DS2yUtTO8+Q5IoquekntjrNZvknK8lScafH/H4yGPdCqyX9ElJM9Oaz4FKh2UX+E/gNJJ+z8u2MY5aPgO8LiIeyXnsuyTNwW9MY51P8j5eAiDpxZLOqJ7rtK/2ncDNI3cUEc+S9Md+RdJb08/5tLTv6osjy1v9nKAsaw5JInqG5NftWuDs9LFPAg8BN6dNLdcBLwKIiAdJOv6vA35J0q+VdSpJ09ETJL9C/6P6QESsBY4CzkiPdyZwVEQ8RfL5/BhJrehpkoEKH8rs9xaSzvOngH8ETkj3B0nH+HSSzvJngMtJamFExGVp+f8E1gPfIxlKDPBPwGfTZqqP552k9AvpwySJ7jGSZL4qr+wEVEcSVv+qHfAHA7dI2gBcCZyW9vdA8sPgwjT2P57g8b9F8hl4jOQcjvxivgA4ID3W99IfIkeRDKL5Ncl7cj5JrafIxSTv6Y/T93tb4igUEauL+ujSvsx3krzfT5PUgm4hGUoOyefi90nO9cb0uPeQfE7z9vcvJJ/Vz5K0GjxK8rn/Xr3x2mgKL1hok5Ckk4D3R4QvhDSbolyDMjOztuQEZWZmbclNfGZm1pZcgzIzs7Y0JS8s22233WLRokWtDsPMzOqwYsWKpyJi1LWEUzJBLVq0iOXLl7c6DDMzq4Ok3KnG3MRnZmZtyQnKzMzakhOUmZm1JScoMzNrS05QZmbWlpygzMysLTlBmZlZW3KCMjMryVMbtvDAE+tbHcak4QRlZlaSr17/MO+78LZWhzFpOEGZmZXkub4BNm4ZaHUYk4YTlJlZSSJgsOIVJOrlBGVmVpJKBM5P9XOCMjMrScU1qHFpSYKS9HZJ90qqSOopKLOXpOsl3ZeWPa3sOM3MGikiGPQisXVrVQ3qHuB44MYaZQaAMyLiAOA1wEckHVBGcGZmzVCJoOIaVN1ash5URKwEkFSrzOPA4+nt9ZJWAguA+8qI0cys0SoBA05QdZsUfVCSFgEHAbfUKHOKpOWSlvf29pYWm5lZvaqpybWo+jQtQUm6TtI9OX/HjnM/s4DvAqdHxLqichGxJCJ6IqJn3rxRKwebmbVcJe1/cj9UfZrWxBcRR050H5KmkSSniyLiiolHZWbWOlFNUJVgWmeLg5kE2raJT0kH1QXAyog4p9XxmJlNVKWS/usaVF1aNcz8OEmrgEOAqyRdnW6fL2lZWuxQ4N3A6yXdmf4tbkW8ZmaNUE1MHihRn1aN4lsKLM3ZvhpYnN6+CSge5mdmNsl4kMT4tG0Tn5nZVJPtg7KxOUGZmZWkmpc8iq8+TlBmZiWp9kFVB0tYbU5QZmYlcQ1qfJygzMxKMtQHNegEVQ8nKDOzkoRrUOPiBGVmVpKKR/GNixOUmVlJhgZJuAZVFycoM7OSDA2ScA2qLk5QZmYl8YW64+MEZWZWknANalycoMzMSuL1oMbHCcrMrCTVipMni62PE5SZWUncBzU+TlBmZiXxVEfj4wRlZlaSwDWo8XCCMjMrSXUWcyeo+jhBmZmVxDNJjI8TlJlZSYavg2ptHJOFE5SZWUk8Wez4OEGZmZVkZBPf3//gPk6/5I5WhtTWulodgJnZ9qJabxpIa1APPrmeZ57ra11Abc41KDOzksSImST6Byvuj6qhJQlK0tsl3SupIqlnjLKdku6Q9IOy4jMza4aRfVD9gzE0u4SN1qoa1D3A8cCNdZQ9DVjZ3HDMzJpv5GSxA4MVDzmvoSUJKiJWRsQDY5WTtBB4C3B+86MyM2uu6oW61Sa+vsHAA/qKtXsf1LnAmcCYrbSSTpG0XNLy3t7epgdmZratqoMkXIOqrWkJStJ1ku7J+Tu2zucfBayJiBX1lI+IJRHRExE98+bNm1DsZmbNMHKYef9gxUtv1NC0YeYRceQEd3EocIykxUA3MEfSdyLiXROPzsysfHmDJDo71MqQ2lrbNvFFxKcjYmFELALeAfzYycnMJrOh5TYyw8zdxFesVcPMj5O0CjgEuErS1en2+ZKWtSImM7Nmi5wmPuenYi2ZSSIilgJLc7avBhbnbL8BuKHpgZmZNdFwDSr5t38wXIOqoW2b+MzMpprhJd+TDOUmvtqcoMzMSjK6BuWpjmpxgjIzK0l2JonBSnKRrqc6KuYEZWZWkuxksf1p1clNfMWcoMzMSpKtQQ0nqFZG1N6coMzMShKZ66AGBrcecm6jOUGZmZUkO5PEUA3KVahCTlBmZiXJ1qD63MQ3JicoM7OSZCeLdRPf2JygzMxKktfE5/xUzAnKzKwk1VxUiaDfNagxOUGZmZUgIoZqSwODwzWoQSeoQk5QZmYlyOah7HVQ4dkkCjlBmZmVINuUl8wkMXzf+SmfE5SZWQkqW9WgGKpBJY85Q+VxgjIzK8HoGlQ2QbUiovbnBGVmVrKBSmWrJj7XoPI5QZmZlSCbhAYrbuKrhxOUmVkJss14lQgGKm7iG4sTlJlZCbauQQX9A27iG4sTlJlZCSKztHslhieLHfmYDXOCMjMrQbB1DWogk6A8m0Q+JygzsxJk+5kGRlyo6ya+fC1JUJLeLuleSRVJPTXK7STpckn3S1op6ZAy4zQza5SR10H1eRTfmFpVg7oHOB64cYxy/xf4n4h4MfByYGWzAzMza4atBklk1oMCT3VUpKsVB42IlQCSCstImgscBpyUPqcP6CshPDOzhssmodEzSThD5WnnPqh9gV7gPyTdIel8STsWFZZ0iqTlkpb39vaWF6WZWR1G1qD6M9dBDfpCqFxNS1CSrpN0T87fsXXuogt4JfD/IuIgYCPwqaLCEbEkInoiomfevHkNeAVmZo2TrSQNDG59HZQrUPma1sQXEUdOcBergFURcUt6/3JqJCgzs3a21SCJcBNfPdq2iS8ingAelfSidNMRwH0tDMnMbJtttWBhxVMd1aNVw8yPk7QKOAS4StLV6fb5kpZliv4FcJGku4FXAJ8vPVgzswbYugYFfZ7qaEzjauKTtDMwH9gEPBKxbRN0RMRSYGnO9tXA4sz9O4HC66TMzCaLyogaVLaJz0u+5xszQaXDvT8CvBOYTjKyrhvYQ9LNwFcj4vqmRmlmNslVk9C0To1q4hv0XHy56qlBXQ58C/jDiPhd9gFJrwLeLen5EXFBE+IzM5sSqjWoro4OBivhJr46jJmgIuINNR5bAaxoaERmZlNQtQbV1alkJomKR/GNpa5BEpK6lE77IGkvSSdIOqi5oZmZTR3VGtS0zo5RM0k4P+UbM0FJ+gCwBvhNevtHwAnAJZI+2eT4zMymhGotqasjqUF5wcKx1dMHdTrwAmA2yWSt+0TEU5J2AG4D/rl54ZmZTQ2RqUFt6h/0VEd1qCdB9UXEM8Azkh6KiKcAIuI5SZ681cysDpVsH1TaxNfVIQYq4Qt1C9SToGam/U0dwPT0ttK/7mYGZ2Y2VcTQKL40QQ0EM7o6GOgb9HVQBepJUE8A5+Tcrt43M7MxVIaug+pI5uKrVJje1cHGvkHXoArUM8z88BLiMDOb0vKa+GZ0dQL9HiRRoJ6ZJI6v9XhEXNG4cMzMpqbshbqVgP6BoHtaMpC64ipUrnqa+I5O/90deC3w4/T+HwH/CzhBmZmNaXiqI4AtA4PMmZl8BTs/5aunie/PACRdAxwQEY+n958HfLOp0ZmZTRHZGhTA5v5qE5+vgyoynuU29qomp9STwN4NjsfMbEqqNuN1pTWozQODzOhKm/icoHKNZ7mNH6XrNl2c3j8RuK7xIZmZTT2VzDBzSIadz0j7oJyf8tWdoCLiVEnHAYelm5ak6zqZmdkYIjPMvKraxOeZJPLVM4pPkZ7ZooUGs2XMzGy06hdkNkHNmzUDcBNfkXr6oK6X9BeStupvkjRd0uslXQi8tznhmZlNDdnroKr22mVm+lhLQmp79TTxvQk4GbhY0r7A74CZJMntGuDciLijaRGamU0BI0fxASzceQfAS74XqWeY+Wbgq8BXJU0DdgM2jVxd18zMilVi6+ugABbu7BpULfUuWNgp6f6I6I+Ix52czMzGp1pL6uwYTlA77TANgEHXoHLVlaAiYhB4YGQ/lJmZ1ae6/FN1kET3tA46VB1y7gSVZzwX6u4M3CvpR5KurP5ty0ElvV3SvZIqknpqlPvLtNw9ki6W5OU9zGxSGh7FlySlPed0DyUoj+LLN54Ldf+qgce9Bzge+HpRAUkLgI+STK+0SdJ/Ae/A0yuZ2SQ0PIovqRfsnk1QlcKnbdfGc6HuTyTtA+wfEdelS753bstBI2IlgKSxinaRLJjYD+wArN6W45mZtVq1Ga9/IMlGe87ppvoV6BpUvrqb+CR9ALic4VrPAuB7TYgJgIh4DPgS8FvgceDZiLimRnynSFouaXlvb2+zwjIz2ybVkXpr1m8BYM+53XRkpj2y0cbTB/UR4FBgHUBE/JJkCY5ckq5L+45G/h1bz8Ek7QwcC+wLzAd2lPSuovIRsSQieiKiZ968eeN4WWZmzVetJU1PJ4jdf/dZdKZVKI/iyzeePqgtEdFXbZaT1MVwv98oEXHkBGM7Evh1RPSmx7uCZD2q70xwv2ZmpavmoJMP3ZdDnr8rxx20gKc2JLUpN/HlG08N6ieSPkPSJ/QG4DLg+80JC0ia9l4jaQclWfEIYGUTj2dm1jTVJDRjWgdve9VCOjo01A/vC3XzjSdBfQroBX4B/DmwDPjsthxU0nGSVgGHAFely3ggab6kZQARcQtJn9ft6TE7gCXbcjwzs1arVpI6MoPDqtfs+jqofONp4vsj4DsR8Y2JHrRoVvSIWA0sztz/G+BvJno8M7NWq9agsmOXh4eZO0HlGU8N6j3AXZJulnS2pKPTgQxmZjaGSl4NqqM6SKIVEbW/8VwH9V5ImuGAE4CvkIyuG08tzMxsu1Rtxste/ukmvtrqTi7pEO8/BF4KPAWcB/y0SXGZmU0pQ31QHdk+KE91VMt4aj/nAg8DXwOuj4hHmhGQmdlUVE1CHVvVoDyKr5a6+6AiYjeShQu7gX+UdKukbzctMjOzKaSahJQZJuGpjmobz1RHc4C9gX2ARcBcwFMcmpnVoWYNylWoXONp4rsp83deRKxqTkhmZlPP8CCJ4QxVXbzQ+SnfeEbxvQxA0qzmhWNmNjVVc1BHzig+N/HlG08T34GS7gDuBe6TtELSgc0Lzcxs6qg242Wvg/JUR7WN50LdJcDHImKfiNgbOANPPWRmVpe8C3WT+74Oqsh4EtSOEXF99U5E3ADs2PCIzMymoKFmvBHrtHZIDLoKlWs8gyR+JemvgOrQ8ncBv2p8SGZmU8/wZLFbb+/okJv4CoynBnUyMA+4AvguUL0uyszMxhCM7oNK7ruJr8iYNShJ3cAHgf1Ilr04IyL6mx2YmdlUUtwHJY/iK1BPDepCoIckOb0ZOLupEZmZTUGVnMlioZqgWhDQJFBPH9QBEfFSAEkXALc2NyQzs6lnaIzEiAQl4UESBeqpQQ0150XEQBNjMTObsvKug4JkNgn3QeWrpwb1cknr0tsCZqb3BUREzGladGZmU0TtPqgWBDQJjJmgIqKzjEDMzKay4VF8W2/vkKc6KjKeYeZmZraNhpbbGFGDkmtQhZygzMxKEBGjak/g66BqcYIyMytBJWJU7Qmg01MdFXKCMjMrQSVG9z+Bm/hqaUmCknS2pPsl3S1pqaSdCsq9SdIDkh6S9KmSwzQza5iI0f1PAB0dbuIr0qoa1LXAgekiiA8Cnx5ZQFIn8BWS2SsOAN4p6YBSozQza5DiPihPdVSkJQkqIq7JXPR7M7Awp9irgYci4lcR0QdcAhxbVoxmZo1UiRh1DRT4Oqha2qEP6mTghznbFwCPZu6vSrflknSKpOWSlvf29jY4RDOzianEqKWggHSqI9egco1nPahxkXQdsGfOQ2dFxH+nZc4CBoCLJnq8iFhCusJvT0+P320zaytFNahOeaqjIk1LUBFxZK3HJZ0EHAUcEfnvzmPAXpn7C9NtZmaTTjJIYvT2DolKpfx4JoNWjeJ7E3AmcExEPFdQ7DZgf0n7SpoOvAO4sqwYzcwaKSLoyBklIU91VKhVfVDnAbOBayXdKelrAJLmS1oGQzOnnwpcDawE/isi7m1RvGZmE5JcB+VBEuPRtCa+WiJiv4Ltq4HFmfvLgGVlxWVm1iyVomHmHa5BFWmHUXxmZlNeUkvKHyThBJXPCcrMrARFF+p6qqNiTlBmZiWIwj4oT3VUxAnKzKwEhX1QbuIr5ARlZlaCStFksV5uo5ATlJlZCSIi/0LdDtwHVcAJysysBLUmi3UfVD4nKDOzEgT5Cxb6Qt1iTlBmZiUomknCUx0Vc4IyMytBpagPSqLiKlQuJygzsxJE0XIbHW7iK+IEZWZWgkqlaLkNN/EVcYIyMytB0Sg+T3VUzAnKzKwEQdGFup7qqIgTlJlZCYomi/VMEsWcoMzMSlB7wUInqDxOUGZmJSgcZt4hnJ/yOUGZmZWgeLJYj+Ir4gRlZlaCWn1Q7oLK5wRlZlaCogULPdVRMScoM7MS1Fyw0FWoXE5QZmYlSAZJ5Ex15Ca+Qk5QZmYlqATkVKDSBQudofK0JEFJOlvS/ZLulrRU0k45ZfaSdL2k+yTdK+m0FoRqZtYQRZPFeqqjYq2qQV0LHBgRLwMeBD6dU2YAOCMiDgBeA3xE0gElxmhm1jARSW1pJE91VKwlCSoiromIgfTuzcDCnDKPR8Tt6e31wEpgQXlRmpk1Tq0l3wedoHK1Qx/UycAPaxWQtAg4CLilRplTJC2XtLy3t7exEZqZTVDxhboexVekq1k7lnQdsGfOQ2dFxH+nZc4iacq7qMZ+ZgHfBU6PiHVF5SJiCbAEoKenx++2mbWViMgfJCFPdVSkaQkqIo6s9bikk4CjgCOioAFW0jSS5HRRRFzR8CDNzEqSTBY7erunOirWtARVi6Q3AWcCr4uI5wrKCLgAWBkR55QZn5lZowUFfVBe8r1Qq/qgzgNmA9dKulPS1wAkzZe0LC1zKPBu4PVpmTslLW5RvGZmE5Is+Z4/1ZEHSeRrSQ0qIvYr2L4aWJzevon869rMzCadoqmOOiUPMy/QDqP4zMymvKLJYj2beTEnKDOzEhQuWOhBEoWcoMzMShAULbeRDDN3M99oTlBmZiUorkEpfbzkgCYBJygzsxIU90El/7qZbzQnKDOzEhQuWNihocdta05QZmYlKFqwsFqrcn4azQnKzKwEyYW6o7e7ia+YE5SZWUmKroMCGPQoiVGcoMzMSjBYye+D0lANqtx4JgMnKDOzEqzb3M/s7mmjtnd2VPugnKFGcoIyM2uyzf2DPNc3yC47Th/1mK+DKuYEZWbWZGs39gGwa26CSv71IInRnKDMzJrs6Q1JgsqrQXVP6wRg45aBUmOaDJygzMyabO3GLQDsOmt0gtpjTjcAT67bUmpMk4ETlJlZkz29sVqDmjHqsT3nJgnqiXWbS41pMnCCMjNrsuEEVVyDWuMENYoTlJlZk63d2EdXh5jTPXoR8zndXXRP6+CJZ52gRnKCMjNrsqc39LHzjtNz5+KTxJ5zut3El8MJysysyZ5+ri93iHnVHnO6WeNBEqM4QZmZNdnTG/ty+5+q9nANKpcTlJlZk42VoPacmyQoT3e0tZYkKElnS7pf0t2SlkraqUbZTkl3SPpBiSGamTXM2g1bxmzi6xuo8Oym/hKjan+tqkFdCxwYES8DHgQ+XaPsacDKUqIyM2ugx5/dxGXLH2Xd5oHca6Cq9pzja6HytCRBRcQ1EVGd1+NmYGFeOUkLgbcA55cVm5lZI0QEp11yJ5+4/G4AdsmZRaJqjzlJ8rrxwV7WrHeSqho9KL98JwOXFjx2LnAmMLu0aMzMaogItgxU2LBlgHWb+nmub5DZ3V08svY5fvdcH9M6O1j+yDP0DQ5y66+f5sAFc7jnsXXsMbu4BrVw5x0A+Pyy+znn2gd52ysXsuecbmZ1dzFrRhezu6cxO709q7uL2TO6mDm9k84ODc2G3iHRoeRfidwh7ZNN0xKUpOuAPXMeOisi/jstcxYwAFyU8/yjgDURsULS4XUc7xTgFIC99957m+P+p2UruXT5o9v8/KpG9HU2qsO0Yd2ujXhNE99Fsp82OjeN6teOBp2dhnz2Jr6Lhu6oEeemUe9TJWLMpTGmd3bQN1jhpQvmcsWHX8vPH17La56/a2H5Ped2c9kHD2HjlgG+e/tjXL5iFVsGKhOKUxpOWpLQiMeGbmce2Xp7trxGb09vvGqfnfnmn716QrEWaVqCiogjaz0u6STgKOCIyP+2ORQ4RtJioBuYI+k7EfGuguMtAZYA9PT0bPNH8aUL57K5f3Bbn76VdvoF06hQtv6Yb+M+GhZLg/bTgB016r1u2CemEa+pQdG00/vdkPcascOMTmantZmZ07pYt7mfBTvNZN7sGWzcMsBL5s/l6Y197Dijk2mdHRz2wnlj7vfgRbsAcPiLdgegb6DCxi0DbNgywPrNyb8btvSzfnNyf1Pf4FCyrEQQEUSw1f3q7cHsV2z+za1+9BUUH9qe/cGw9y47jPnatpVaMaxR0puAc4DXRURvHeUPBz4eEUfVs/+enp5Yvnz5hGI0M7NySFoRET0jt7dqFN95JP1K10q6U9LXACTNl7SsRTGZmVkbackgiYjYr2D7amBxzvYbgBuaG5WZmbUTzyRhZmZtyQnKzMzakhOUmZm1JScoMzNrS05QZmbWlpygzMysLbXkQt1mk9QL/GYCu9gNeKpB4Ux2PhcJn4eEz8Mwn4tEI87DPhExarqNKZmgJkrS8ryrmrdHPhcJn4eEz8Mwn4tEM8+Dm/jMzKwtOUGZmVlbcoLKt6TVAbQRn4uEz0PC52GYz0WiaefBfVBmZtaWXIMyM7O25ARlZmZtabtOUJLeJOkBSQ9J+lTO4zMkXZo+foukRS0Is+nqOA+HSbpd0oCkE1oRY1nqOBcfk3SfpLsl/UjSPq2Is9nqOA8flPSLdD23myQd0Io4m22s85Ap9zZJIWnKDjuv4zNxkqTe9DNxp6T3T/igMbRU8Pb1B3QCDwPPB6YDdwEHjCjzYeBr6e13AJe2Ou4WnYdFwMuAbwEntDrmFp+LPwJ2SG9/aDv+TMzJ3D4G+J9Wx92K85CWmw3cCNwM9LQ67hZ+Jk4CzmvkcbfnGtSrgYci4lcR0QdcAhw7osyxwIXp7cuBIySpxBjLMOZ5iIhHIuJuoNKKAEtUz7m4PiKeS+/eDCwsOcYy1HMe1mXu7ghMxdFW9XxHAPw98M/A5jKDK1m956KhtucEtQB4NHN/Vbott0xEDADPAruWEl156jkP24vxnov3AT9sakStUdd5kPQRSQ8DXwQ+WlJsZRrzPEh6JbBXRFxVZmAtUO//jbelzd+XS9progfdnhOU2TaT9C6gBzi71bG0SkR8JSJeAHwS+Gyr4ymbpA7gHOCMVsfSJr4PLIqIlwHXMtz6tM225wT1GJDN8AvTbbllJHUBc4G1pURXnnrOw/airnMh6UjgLOCYiNhSUmxlGu9n4hLgrc0MqEXGOg+zgQOBGyQ9ArwGuHKKDpQY8zMREWsz/x/OB1410YNuzwnqNmB/SftKmk4yCOLKEWWuBN6b3j4B+HGkvYFTSD3nYXsx5rmQdBDwdZLktKYFMZahnvOwf+buW4BflhhfWWqeh4h4NiJ2i4hFEbGIpE/ymIhY3ppwm6qez8TzMnePAVZO9KBdE93BZBURA5JOBa4mGaHy7xFxr6S/A5ZHxJXABcC3JT0EPE3ypkwp9ZwHSQcDS4GdgaMl/W1EvKSFYTdFnZ+Js4FZwGXpeJnfRsQxLQu6Ceo8D6emNcl+4BmGf8hNGXWeh+1Cnefio5KOAQZIvi9PmuhxPdWRmZm1pe25ic/MzNqYE5SZmbUlJygzM2tLTlBmZtaWnKDMzKwtbbfDzM0aSdIg8IvMprdGxCMtCsdsSvAwc7MGkLQhImYVPCaS/2tTfbJds4ZyE59ZE0halK6d8y3gHmAvSZ+QdFs6mebfZsqeJenBdF2liyV9PN1+Q3XaHEm7pdPpIKlT0tmZff15uv3w9DmXS7pf0kXV2fclHSzpfyXdJelWSbMl3SjpFZk4bpL08rLOkdlY3MRn1hgzJd2Z3v418JfA/sB7I+JmSW9M778aEMmcbYcBG0lmKHkFyf/H24EVYxzrfcCzEXGwpBnAzyRdkz52EPASYDXwM+BQSbcClwInRsRtkuYAm0hmSjkJOF3SC4HuiLhrYqfBrHGcoMwaY1NEvKJ6R8nqy7+JiJvTTW9M/+5I788iSVizgaXVNaYk1TN9zhuBl2l4deO56b76gFsjYlW6rztJFpt8Fng8Im6D4bWcJF0G/JWkTwAnA98c52s2ayonKLPm2Zi5LeCfIuLr2QKSTq/x/AGGm+G7R+zrLyLi6hH7OhzIzq4+SI3/4xHxnKRrSRae+2MaMPu0WSO5D8qsHFcDJ0uaBSBpgaTdSZYKf6ukmZJmA0dnnvMIw0njhBH7+pCkaem+XihpxxrHfgB4XjrpL2n/UzVxnQ/8K3BbRDwzoVdo1mCuQZmVICKukfR7wM/TcQsbgHdFxO2SLgXuAtaQLGtQ9SXgvySdAmRXbD2fpOnu9nQQRC811mOKiD5JJwL/JmkmSf/TkcCGiFghaR3wH415pWaN42HmZm1E0udIEseXSjrefOAG4MUeBm/txk18ZtspSe8BbgHOcnKyduQalJmZtSXXoMzMrC05QZmZWVtygjIzs7bkBGVmZm3JCcrMzNrS/wcTj94ZpcO4ngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.linalg import toeplitz\n",
    "import math\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "N=256\n",
    "n=np.arange(256)\n",
    "x=(np.exp(2*np.pi*0.4*n*1j)+np.exp(2*np.pi*0.425*n*1j)+\n",
    "   0.2/math.sqrt(2)*np.random.normal(size=n.shape)+np.random.normal(size=n.shape)*1.0j)\n",
    "\n",
    "# 自相关函数\n",
    "def xcorr(data):\n",
    "    length =len(data)\n",
    "    R=[]\n",
    "    for m in range (0,length):\n",
    "        sums=0.0\n",
    "        for nn in range(0,length-m):\n",
    "            sums=sums+data[nn]*data[nn+m]\n",
    "        R.append(sums/length)\n",
    "    return R\n",
    "# 求自相关\n",
    "R=xcorr(x)\n",
    "# 自相关矩阵\n",
    "Rx =toeplitz(R)\n",
    "w,v = np.linalg.eig(Rx)\n",
    "v=np.mat(v)\n",
    "vh=v.H\n",
    "P=[]\n",
    "for index in range(0, N+1):\n",
    "    ew = []\n",
    "    w = index / N * np.pi\n",
    "    for k in range(0, N):\n",
    "        item = complex(np.cos(k * w), np.sin(k * w))\n",
    "        ew.append(item)\n",
    "    ew = np.mat(ew)\n",
    "    ew = ew.T\n",
    "    eHw = ew.H\n",
    "    sums = 0\n",
    "    sums = np.mat(sums)\n",
    "    for j in range(3, N + 1):\n",
    "        sums = sums + v[:, j - 1] * vh[j - 1,:]\n",
    "    fenzi = eHw * sums * ew\n",
    "    B=1/fenzi[0,0]\n",
    "    B=math.log10(B)\n",
    "    P.append(B)\n",
    "\n",
    "f = np.linspace(0, N, N + 1) / (2 * N)\n",
    "plt.plot(f,P)\n",
    "plt.title('Pseudospectrum Estimate via MUSIC')\n",
    "plt.xlabel(\"Frequency\")\n",
    "plt.ylabel(\"Power(dB)\")\n",
    "plt.tight_layout() # 自动布局\n",
    "plt.show"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
